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qq ■ When hadrons scatter at high energies, strong color fields, whose dynamics is de- 

scribed by quantum chromodynamics (QCD), are generated at the interaction point. 
If one represents these fields in terms of partons (quarks and gluons), the average 

■ number densities of the latter saturate at ultrahigh energies. At that point, non- 
linear effects become predominant in the dynamical equations. The hadronic states 

Qh' that one gets in this regime of QCD are generically called "color glass condensates" . 

Our understanding of scattering in QCD has benefited from recent progress in 
statistical and mathematical physics. The evolution of hadronic scattering ampli- 
tudes at fixed impact parameter in the regime where nonlinear parton saturation 
effects become sizable was shown to be similar to the time evolution of a system of 
\ classical particles undergoing reaction-diffusion processes. The dynamics of such a 

system is essentially governed by equations in the universality class of the stochastic 
Fisher-Kolmogorov-Petrovsky-Piscounov equation, which is a stochastic nonlinear 
partial differential equation. Realizations of that kind of equations (that is, "events" 
in a particle physics language) have the form of noisy traveling waves. Universal 
properties of the latter can be taken over to scattering amplitudes in QCD. 

This review provides an introduction to the basic methods of statistical physics 
useful in QCD, and summarizes the correspondence between these two fields and 

H ' its theoretical and phenomenological implications. 
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1 Introduction to high energy scattering in QCD 



The study of quantum chromodynamics in the high-energy regime has taken 
a new soar in the last 15 years with the wealth of experimental data that 
have been collected, first at the electron-proton collider DESY-HERA, and 
then at the heavy-ion collider RHIC. More energy in the collision enables the 
production of objects of higher mass in the final state, and thus the discovery 
of new particles. But higher energies make it also possible to observe more 
quantum fluctuations of the incoming objects, that is to say, to study more 
deeply the structure of the vacuum. 

The well-established microscopic theory which describes the interactions of 
hadronic objects is quantum chromodynamics (QCD). (For a comprehensive 
textbook, see Ref. [1]). There are not many known analytical approaches to 
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QCD, except perturbative expansions of observables in powers of the strong 
coupling constant a s which, thanks to asymptotic freedom, is justified for 
carefully chosen observables in special kinematical regimes. But fixed-order 
calculations in QCD are known to usually have a very limited range of ap- 
plicability. This is because in the evaluation of Feynman graphs, the coupling 
constant always comes with "infrared" and "collinear" logarithms that are 
related to the phase space that is available to the reaction, that is to say, to 
kinematics. Resumming part of these logarithms is mandatory. All of them is 
too difficult. The question is to carefully select the dominant ones, and this is 
not at all easy. 

At the HERA collider, electrons or positrons scattered off protons at the 
center-of-mass energy y/s, exchanging a photon of virtuality Q. Through the 
scattering, one could probe partonic fluctuations of the proton (made of quarks 
and gluons) of transverse momenta k ~ Q, and longitudinal momentum frac- 
tions x ~ Q 2 /(Q 2 + s). 

For a long time, the dominant paradigm had been that the collinear logarithms 
logQ 2 , that become large when Q 2 is large compared to the QCD confinement 
scale A 2 , were the most important ones. As a matter of fact, searches for 
new particles or for exotic physics require to scrutinize matter at very small 
distances, and hence very large Q 2 have to be considered. Perturbative series 
of powers of a s log Q 2 have to be fully resummed. The equation that performs 
this resummation is the celebrated Dokshitzer-Gribov-Lipatov-Altarelli-Parisi 
(DGLAP) equation [2,3,4]. 

However, once HERA had revealed its ability to get extremely good statistics 
in a regime in which Q 2 is moderate (from 1 to 100 GeV 2 ) and x very small 
(down to 10~ 5 ) it became clear that infrared logarithms (log 1/x) could show up 
and even dominate the measured observables. The resummation of the series 
of infrared logs is performed by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) 
equation [5,6,7]. The series X)(a: s 1°§ ^/ x ) k (with appropriate coefficients) is the 
leading order (LO), while the series ^a s (a s \ogl/x) k is the next-to-leading 
order (NLO), which has also been computed [8,9]. The BFKL equation is a 
linear integro-differential equation. 

At ultrahigh energy, the bare BFKL equation seems to violate the Froissart 
bound, that states that total hadronic cross sections cannot rise faster than 
(log 2 s)/m 2 . The latter is a consequence of the unitarity of the probability of 
scattering. The BFKL equation predicts a power rise with the energy of the 
form s £ , where e is positive and quite large (0.3 to 0.5 according to the effective 
value of a s that is chosen). The point at which the BFKL equation breaks down 
depends on the value of the typical transverse momentum which characterizes 
the observable (It is the photon virtuality Q in the case of deep-inelastic 
scattering). One may define the energy-dependent saturation scale Q s (x) in 
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such a way that the BFKL equation holds for Q > Q s (x). For Q ~ Q s (x), 
the probability for scattering to take place is of order 1, and for Q < Q s (x), 
it would be larger than 1 if one trusted the BFKL equation. The saturation 
scale is a central observable, which we shall keep discussing in this review: It 
signs the point at which the linear (BFKL) formalism has to be corrected for 
nonlinear effects. The regime in which nonlinearities manifest themselves is a 
regime of strong color fields, sometimes called the color glass condensate (For 
the etymology of this term, see e.g. the lectures of Ref. [10]; for a review, see 
Ref. [11]). 

The fact that unitarity is violated is not only due to the lack of a hadronic 
scale in the BFKL equation, which is a perturbative equation; Introducing 
confinement in the form of a cutoff would not help this particular problem. It 
simply means that still higher orders are needed. The NLO corrections to the 
BFKL kernel indeed correct this behavior in such a way that the description 
of the HERA data in the small- x regime is possible by the BFKL equation. 
However, these corrections are not enough to tame the power-like growth 
of cross sections as predicted by the LO BFKL equation. It seems that a 
resummation of contributions of arbitrary order would be needed. 

New equations were proposed well before the advent of colliders able to reach 
this regime. Gribov-Levin-Ryskin wrote down a model for the evolution of 
the hadronic scattering cross sections in the early 80's [12,13], and Mueller 
and Qiu derived a similar equation from QCD a bit later [14]. These equa- 
tions are integral evolution equations with a nonlinear term, which basically 
takes into account parton saturation effects, that is to say, recombination or 
rescattering. The latter cannot be described in a linear framework such as the 
BFKL formalism. Subsequently, more involved QCD evolution equations were 
derived from different points of view. In the 90 's, McLerran and Venugopalan 
[15,16,17] proposed a first model, mainly designed to approach heavy-ion colli- 
sions. Subsequently, Balitsky [18], Jalilian-Marian, Iancu, McLerran, Weigert, 
Leonidov and Kovner (B-JIMWLK) [19,20,21,22,23] worked out QCD correc- 
tions to this model, and got equations that reduce to the BFKL equation 
in the appropriate limit. Technically, these equations actually have the form 
of an infinite hierarchy of coupled integro-differential equations (in Balitsky's 
formulation [18]), of a functional renormalizaton group equation, or alterna- 
tively, of a Langevin equation (in Weigert's formulation [23]). A much simpler 
equation was derived in 1996 by Balitsky [18] and rederived by Kovchegov in 
1999 [24,25] in a very elegant way within a different formalism. The obtained 
equation is called the Balitsky-Kovchegov equation (BK). The latter deriva- 
tion was based on Mueller's color dipole model [26], which proves particularly 
suited to represent QCD in the high energy limit. 

The exciting feature of this kinematical regime of hadronic interaction from a 
theoretical point of view is that the color fields are strong, although, at suf- 
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ficiently high energies, the QCD coupling is weak, authorizing a perturbative 
approach, and thus analytical calculations. In such strong field regime, nonlin- 
ear effects become crucial. But the conditions of applicability of the different 
equations that had been found had never been quite clear. Anyway, these 
equations are extremely difficult to solve, which had probably been the main 
obstacle to more rapid theoretical developments in the field until recently. 

Furthermore, for a long time, the phenomenological need for such a sophisti- 
cated formalism was not obvious, since linear evolution equations such as the 
DGLAP equation were able to account for almost all data. But Golec-Biernat 
and Wusthoff showed that unitarization effects may have already been seen at 
HERA [27,28]. Their model predicted, in particular, that the virtual photon- 
proton cross section should only depend on one single variable r, made of a 
combination of the transverse momentum scale (fixed by the virtuality of the 
photon Q) and x. This phenomenon was called "geometric scaling" [29]. It 
was found in the HERA data (see Fig. 1): This is maybe one of the most 
spectacular experimental result from HERA in the small-x regime. 

This observation has triggered many phenomenological and theoretical works. 
Soon after its discovery in the data, geometric scaling was shown to be a so- 
lution of the Balitsky-Kovchegov (BK) equation, essentially numerically, with 
some analytical arguments (see e.g. [31,32,33,34]). The energy dependence of 
the saturation scale was eventually precisely computed by Mueller and Tri- 
antafyllopoulos [35]. Later, it was shown that the BK equation is actually in 
the universality class of the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) 
equation [36,37], and geometric scaling was found to be implied by the fact 
that the latter equation admits traveling wave solutions [38]. 

A first step beyond the BK equation, in the direction of a full solution to high 
energy QCD, was taken by Mueller and Shoshi in 2004 [39]. Actually, they 
did not solve the B-JIMWLK equations, but instead, they solved the linear 
BFKL equation with two absorptive boundary conditions, which they argued 
to be appropriate to represent the expected nonlinearities. Geometric scaling 
violations were found from their calculation, which should show up at any 
energies. 

Subsequently, it was shown that high-energy QCD at fixed coupling is actually 
in the universality class of reaction- diffusion processes, studied in statistical 
physics, whose dynamics may be encoded in equations similar to the stochastic 
FKPP equation [40]. The Mueller- Shoshi solution was shown to be consistent 
with solutions to the latter equation. So high-energy QCD seems to be in 
correspondence with disordered systems studied in statistical physics. This 
correspondence has provided a new understanding of QCD in the high-energy 
regime, and it has proven very useful to find more features of high-energy scat- 
tering. The obtained results go beyond a solution to the B-JIMWLK equation, 
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Fig. 1. [From Ref. [30]] Photon-proton total cross section from the most recent 
set of deep-inelastic scattering data in the low-x regime plotted as a function of a 
single scaling variable r = Q 2 /Q 2 (x), where Q is the virtuality of the photon and 
Q 2 (x) ~ A 2 x -0 ' 3 is the so-called saturation scale. Although the cross section is a 
priori a function of two variables, all data fall on the same curve. This phenomenon 
is called geometric scaling [29]. 

which in fact, thanks to the new picture, is seen to be incomplete. 

Scope 

The goal of this review article is to summarize the main ideas behind this 
conjectured correspondence between scattering at high-energy in QCD and 
some processes studied in statistical physics, as well as to introduce the QCD 
reader to the useful technical tools borrowed from statistical physics. We also 
feel that there is a cultural gap to be filled between statistical physics and 
particle physics. Indeed, statistical physicists are used to build simple toy 
models which contain the interesting physics, and whose main properties are 
likely to be independent of the details of the model, i.e. universal. In QCD, 
since the theory is well-established, we are often reluctant to give up some of 
its features to work out exact results in a toy model. One of our aims is to 
convince the reader that such a way of thinking is efficient in the case of high- 
energy QCD, by exhibiting results for QCD scattering amplitudes, obtained 
by looking for the universality class of the considered process, and that are 
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believed to be exact. 

Over the last few years, several hundreds of papers have appeared related to 
this subject, mainly issued from a very active though restricted community. 
Obviously, we cannot give a complete account of this abundant literature. As 
a matter of fact, some important recent developments had to be left out, for 
which we shall only provide references for the interested reader who might 
want to deepen his study in these directions. Concerning the correspondence 
itself, we do not attempt to establish a definite stochastic nonlinear evolution 
equation for QCD amplitudes, for to our judgement, this research line is not 
mature enough yet: A better understanding of the very saturation mecha- 
nism at work in QCD is definitely needed before one may come to this issue. 
Furthermore, it is not clear to us that a stochastic formulation would be a 
technical progress, since there are not many known methods to handle com- 
plicated stochastic equations. We feel that the same is true for the search for 
effective actions that would include so-called Pomeron loops. We also do not 
address the developments based on the boost-invariance symmetry that scat- 
tering amplitudes should have: This would drive us too far off the main focus 
of this review. As for more phenomenological aspects, we only discuss the basic 
features of total cross sections without attempting to address other observ- 
ables such as diffraction. We do also not address the issue of next-to-leading 
effects such as the running of the QCD coupling. This discussion, though 
crucial if one wants to make predictions for actual colliders, would probably 
only be technical in its nature: There is no conceptual difference between the 
fixed coupling and the running coupling case. Here, only basic phenomenolog- 
ical facts brought about by this new understanding of high-energy QCD are 
addressed, namely geometric scaling and diffusive scaling. 

Outline 

The outline goes as follows. The next section is devoted to describing scattering 
in QCD from a s-channel point of view, relying essentially on the parton 
model or rather on an interpretation useful in the high-energy limit, the color 
dipole model. Once this picture is introduced, it is not difficult to understand 
the correspondence with reaction-diffusion processes occuring in one spatial 
dimension, whose dynamics is captured by equations in the university class 
of the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation. We then 
explain how traveling waves appear in this context. In Sec. 3, we study in 
greater detail a toy model for which many technics (field theory, statistical 
methods) may be worked out completely. This model however ignores spatial 
dimensions, and thus, does not account for traveling waves. We summarize the 
state-of-the-art research on equations in the universality class of the FKPP 
equation in Sec. 4. Finally, we come back to QCD, discussing the relevance 
of one-dimensional-like models in the FKPP class, and showing how noisy 
traveling waves may show up in the actual data. 
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2 Hadronic interactions in a s-channel picture and analogy with 
reaction-diffusion processes 

In this section, we shall introduce the physical picture of high-energy scattering 
in the parton model. In particular, the color dipole model [26] is described since 
it is particularly suited to address high-energy scattering, especially close to 
the regime in which nonlinear effects are expected to play a significant role. 
In a second part, we shall argue that high-energy scattering is a peculiar 
react ion- diffusion process. 

2.1 Parton model and dipoles 

2.1.1 General picture 

For definiteness, let us consider the scattering of a hadronic probe off a given 
target, in the restframe of the probe and at a fixed impact parameter, that 
is to say, at a fixed distance between the probe and the center of the target 
in the two-dimensional plane transverse to the collision axis. In the parton 
model, the target interacts through one of its quantum fluctuations, made of 
a high occupancy Fock state if the energy of the reaction is sufficiently high 
(see Fig. 2a). As will be understood below, the probe effectively "counts" the 
partons in the current Fock state of the target whose transverse momenta k 
(or sizes r ~ are of the order of the one that characterizes the probe: 

Roughly speaking, the amplitude for the scattering off this particular partonic 
configuration is proportional to the number of such partons. 

The observable that is maybe the most sensitive to quantum fluctuations of 
a hadron is the cross section for the interaction of a virtual photon with a 
hadronic target such as a proton or a nucleus. The virtual photon is emitted 
by an electron (or a positron). What is interesting with this process, called 
"deep-inelastic scattering", is that the kinematics of the photon is fully con- 
trolled by the measurement of the scattered electron. The photon can be 
considered a hadronic object since it interacts through its fluctuations into a 
quark- ant iquark state. The latter form color dipoles since although both the 
quark and the antiquark carry color charge, the overall object is color neu- 
tral due to the color neutrality of the photon. The probability distribution 
of these fluctuations may be computed in quantum electrodynamics (QED). 
Subsequently, the dipole interacts with the target by exchanging gluons. The 
dipole-target cross section factorizes at high energy. One typical event is de- 
picted in Fig. 2a. 

Dipole models [41,42] have become more and more popular among phenome- 
nologists since knowing the dipole cross section enables one to compute dif- 
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ferent kinds of observables. Like parton densities, the latter is a universal 
quantity, that may be extracted from one process and used to predict other 
observables. Different phenomenological models may be tried for the dipole 
cross section. QCD evolution equations may even be derived, as we shall dis- 
cover below. An accurate recent study of the foundations of dipole models 
may be found in Ref. [43,44]. 

In QCD, the state of a hadronic object, encoded in a set of wave functions, is 
built up from successive splittings of partons starting from the valence struc- 
ture. This is visible in the example of Fig. 2a: The quark and the antiquark 
that build up, in this example, the target in its asymptotic state each emit 
a gluon, which themselves emit, later on in the evolution, other gluons. As 
one increases the rapidity y by boosting the target, the opening of the phase 
space for parton splittings makes the probability for high occupation numbers 
larger. Indeed, the probability to find a gluon that carries a fraction z (up 
to dz) of the momentum of its parent parton (which may be a quark or a 
gluon) is of order a s N c dz/z for small z. There is a logarithmic singularity in 
z, meaning that emissions of very soft gluons (small z) are favored if they are 
allowed by the kinematics. The splitting probability is of order 1 when the 
total rapidity of the scattering y = log \jx is increased by roughly 1/a, where 
the convenient notation a = a s N c / 7r has been introduced. Only splittings of a 
quark or of a gluon into a gluon exhibit the \ j z singularity. Therefore, at large 
rapidities, gluons eventually dominate the partonic content of the hadrons. 

The parton model in its basic form, where the fundamental objects of the 
theory (quarks and gluons) are directly considered, is not so easy to handle in 
the high-energy regime. One may considerably simplify the problem by going 
to the limit of large number of colors (N c ), in which a gluon may be seen as 
a zero-size quark- antiquark pair. Then, color-neutral objects become collec- 
tions of color dipoles, whose endpoints consist in "half gluons" (see Fig. 2b). 
There is only one type of objects in the theory, dipoles, which simplifies very 
much the picture. Furthermore, going to transverse coordinate space (instead 
of momentum space, usually used in the DGLAP formalism) by trading the 
transverse momenta of the gluons for the sizes of the dipoles (through an 
appropriate Fourier transform) brings another considerable simplification. In- 
deed, the splittings that contribute to amplitudes in the high-energy limit are 
the soft ones, for which the emitted gluons take only a small fraction of the 
momentum of their parent (the latter being very large). Therefore, the posi- 
tions of the gluons (and thus of the edges of the dipoles) in the plane transverse 
to the collision axis are not modified by subsequent evolution once the glu- 
ons have been created. Thus, the evolution of each dipole proceeds through 
completely independent splittings to new dipoles. 

We will now see how this picture translates into a QCD evolution equation 
for scattering amplitudes, first in the regime in which there are no nonlinear 
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(a) 



(b) 



Fig. 2. (a) The scattering of a virtual photon probe off a particular fluctuation of 
an evolved target made of a quark and an antiquark in its bare state. The photon 
necessarily goes through a quark- antiquark pair at high enough energies, when the 
target is dominated by dense gluonic states. (What is represented in this figure is 
actually the inelastic amplitude, which is a cut of the total cross section or of the 
forward elastic amplitude), (b) In the dipole model, the probe and the target may 
be represented by a set of color dipoles, and the interaction proceeds through gluon 
exchanges. The curly vertical line represents the lowest-order interaction between 
pairs of dipoles, that is to say, the exchange of a gluon (or a two-gluon exchange if 
one is speaking of the forward elastic amplitude). 

effects. In a second step, we will try and understand how to incorporate the 
latter. 



2.1.2 BFKL equation from the dipole model 

The building up of the states of each hadron is specified by providing the 
splitting rate of a dipole whose endpoints have transverse coordinates (xq, xi) 
into two dipoles (x ,x 2 ) and (xi,x 2 ) as the result of a gluon emission at 
position X2- It is computed in perturbative QCD and reads [26] 

dP 

((x , Xi) -> (x Q , x 2 ), (x 2 , Xi)) 





x 




X\ 


2 


d 2 x 2 


\XQ - 


x 2 
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X\ 


- x 2 


2 2tt 



d(ay) 

Dipole splittings are independent. After some rapidity evolution starting from 
a primordial dipole, one gets a chain of dipoles such as the one depicted in 
Fig. 3. 

The elementary scattering amplitude for one projectile dipole (xq,xx) off a 
target dipole (z , z\) is independent of the rapidity and reads [26] 

n(* , si), (.o, zi)) = ^ log 2 \ X ° ~ Zll 2 Xl ~ Z °l (2) 
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Fig. 3. Schematic picture of a realization of the dipole evolution after the first two 
steps of the evolution ((a) and (b)), and after some larger rapidity evolution (c). 
In the first step (a), the initial dipole {xq,x\) (denoted by a dashed line) splits 
to the new dipoles (xq,X2) and (x2,x±) (full lines). The points represent the edges 
of each dipole, that is to say, the position of the gluons. In the next step (b), the 
dipole (x2,x\) itself splits in two new dipoles. The splitting process proceeds (c) 
until the maximum rapidity is reached. Many very small dipoles are produced in the 
vicinity of each of these endpoints, due to the infrared singularity visible in Eq. (1) 
(Only a fraction of them is represented). The zones 1 and 2 in (c), separated by the 
transverse distance Ab, would evolve quasi-independently after the stage depicted 
in this figure. 
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If the target is an evolved state at rapidity y, then it consists instead in a 
distribution n(y, (z , z±)) of dipoles. The (forward elastic) scattering amplitude 
A(y, (xo,Xi)) is just given by the convolution of n and T cl , namely 

A(y, (x , xi)) = J < ^ < ^-T c \(x , x x ), (z , (z , z ± )). (3) 



Let us examine the properties of T cl . To this aim, it is useful to decompose the 
coordinates of the dipoles in their size vector r a = xq — x\ (resp. = zq — z\) 
and impact parameter b a = X °+ Xl (resp. bb = Z "+ Zl ), In the limit in which the 
relative impact parameters of the dipoles b = b a — bb is very large compared 
to their sizes, we get the simplified expression 

and thus the scattering amplitude decays fast as a function of the relative 
impact parameter. If instead the relative impact parameter is small (of the 
order of the size of the smallest dipole), we get 

T c \r a ,r b ,b) ~ a*jf, (5) 

where r< = min(|r a |, |r&|), r> = max(|r a |, |r(,|), and the integration over the 
angles has been performed. 

Equation (4) means that the dipole interaction is local in impact parameter: 
It vanishes as soon as the relative distance of the dipoles is a few steps in units 
of their size. Eq. (5) shows that only dipoles whose sizes are of the same order 
of magnitude interact. These properties are natural in quantum mechanics. 
Thus the amplitude A in Eq. (3) effectively "counts" the dipoles of size of the 
order of \x i\ at the impact parameter X °+ Xl (up to |x i|), with a weight factor 
a 2 s . 

An evolution equation for the amplitude A with the rapidity of the scatter- 
ing can be established. It is enough to know how the dipole density in the 
target evolves when rapidity is increased, since all the rapidity dependence is 
contained in n in the factorization (3), and such an equation may easily be 
worked out with the help of the splitting rate distribution (1). It reads [26] 

dn(y, (x ,xi)) 
d(ay) 



I — it? (x ,x 2 )) +n{y, (x 2 ,x 1 )) 

J ZTT |£02| F12I 

-n(y, (x ,xi))], (6) 



where x a b = x a —x b . The very same equation holds for A. The elementary scat- 
tering amplitude T el only appears in the initial condition at y — 0, which is not 
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shown in Eq. (6). In a nutshell, the integral kernel encodes the branching dif- 
fusion of the dipoles. The total number of dipoles at a given impact parameter 
grows exponentially, and their sizes diffuse. The appropriate variable in which 
diffusion takes place is log(l/|a;oi| 2 ). (This is due to the collinear singularities 
in Eq. (1).) This equation is nothing but the BFKL equation. A complete so- 
lution to this equation, including the impact-parameter dependence, is known 
[45]. 

An important property of the amplitude A is that it is boost-invariant. This 
property is preserved in the BFKL formulation. We could have put the evolu- 
tion in the projectile instead of the target, or shared it between the projectile 
and the target: The result for the scattering amplitude would have been the 
same. In a frame in which the target carries y' units of rapidity and the pro- 
jectile y — y', the amplitude A reads 

A{y,(xo,xi)) = J ^^^^n^^y-y'^zo^^x,)) 

n xT e \(z ,z 1 ),(ziz[))n^(y>,(z> ,z' 1 )). (7) 

n projcctiic^ _ yi^ ZQjZl ^xo,xi)) is the density of dipoles (zo,zi) found in a 
dipole of initial size (xq,Xi) after evolution over y — y' steps in rapidity. If 
y' = y, one recovers Eq. (3). If y' = 0, then all the evolution is in the projectile 
instead. 

The amplitude A is related to an interaction probability, and thus, it must be 
bounded: In appropriate normalizations, A has to range between and 1. But 
as stated above, the BFKL equation predicts an exponential rise of A with 
the rapidity for any dipole size, which at large rapidities eventually violates 
unitarity. Hence the BFKL equation does not provide a complete account of 
high-energy scattering in QCD. 

2.1.3 Unitarity and the Balitsky-Kovchegov equation 

It is clear that one important ingredient that has been left out in the derivation 
of the BFKL equation is the possibility of multiple scatterings between the 
probe and the target. Several among the n dipoles in Eq. (7) may actually 
interact with the dipoles in the other hadron simultaneously. The only reason 
why such interactions may not take place is that T el ~ a 2 (see Eq. (2)), and 
thus the probability for two simultaneous scatterings is of order a*, which 
is suppressed. But this argument holds only as long as the dipole number 
densities are of order I. If n ~ l/a 2 s (which is also the point above which the 
unitarity of A is no longer preserved in the BFKL approach), then it is clear 
that multiple scatterings should occur. 

In order to try and implement these multiple scatterings, we introduce the 
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S(y,(x ,x 2 )) 



S(y,(x 2 ,X!)) 



Fig. 4. Derivation of the Balitsky-Kovchegov equation. 



probability that there be no scattering between a dipole (xq, x\) and a given re- 
alization of the target with total rapidity y, that we shall denote by S(y, (xq, xi)) 
Let us start with a system in which the evolution is fully contained in the tar- 
get. We increase the total rapidity by boosting the projectile (initially at rest) 
by a small amount dy. Then there are two cases to distinguish, depending on 
whether the dipole (xq,Xi) splits in the rapidity interval dy. In case it splits 
into two dipoles (x ,x 2 ) and (x 2 ,Xi), the probability that the projectile does 
not interact is just the product of the probabilities that each of these new 
dipoles do not interact. This is because once created, dipoles are supposed to 
be independent. In summary: 



S(y + dy, (x ,xi)) = < 



S(y, (x ,xi)) 

with proba 1 - ady J X2 ^gy 
S(y, (x ,x 2 )) x S(y, (x 2 ,x 1 )) 



%01 — ► Xq2, Xl2) 



with proba ady-^g^ 



v Xqi — > Xq2, Xl2) 



Taking the average over the realizations of the target and the limit dy — > 0, 
we get 



— (S(y,(x ,x 1 ))) = a / — ^-[(S(y, (x ,x 2 ))S{y, (x 2 ,x 1 ))) 

ay J Z7T Xq 2 x 21 

^{S(y,(x ,x 1 )))} (9) 

(See Fig. 4 for a graphical representation.) We see that this equation is not 
closed: An evolution equation for the correlator (S(y, (x , x 2 ))S(y, (x 2 , Xi))) is 
required. However, we may assume that these correlators factorize 

(S(y, (x ,x 2 ))S(y, (x 2 ,x 1 ))) = (S(y, (x , x 2 )))(S(y, (x 2 ,xi))). (10) 

This assumption is justified if the dipoles scatter off independent targets, for 
example, off the nucleons of a very large nucleus. Writing A = 1 — (S), we get 
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Fig. 5. Picture of the BK evolution. All the QCD evolution is put in the probe, which 
carries the total rapidity. It develops a high occupancy state of dipoles, which scatter 
independently off the target. 

the following closed equation for A: 

■E-My, {xq, a?i)) = a / — j^-[A(y, (x , x 2 )) + A(y, (x 2 , x x )) 

- A(y, (x , xi)) - A(y, (x , x 2 ))A(y, (x 2 , xi))], (11) 

which is the Balitsky-Kovchegov (BK) equation [24,25]. Note that if one ne- 
glects the nonlinear term, one gets back the BFKL equation (6) (written for 
A instead of n). A graphical representation of this equation is given in Fig. 5. 

It is not difficult to see analytically that the BK equation preserves the uni- 
tarity of A: When A becomes of the order of 1, then the nonlinear term gets 
comparable to the linear terms in magnitude, and slows down the evolution 
of A with y, which otherwise would be exponential. 

Let us go back to Eqs. (8), (9) and instead of assuming the factorization of 
the correlators (10), work out an equation for the two-point correlator (SS). 
From the same calculation as before, we get 

d f d?X3 x'q 2 

t-{So 2 S 2 ii) — a — g — ((S 0S S S2 S 2 'i) — (S 02 S 2 >i)) 

oy J 2vr x z 03 xi 2 

+ a I d l X3 2 \ {(SbsSuSm) - (S 02 S 2 n)) , (12) 
J Zir x 13 x 32 > 

where we have introduced the notation S a b = S(y, (x a , xj)). (See Fig. 6a for 
the corresponding graphical representation.) 

This equation calls for a new equation for the 3-point correlators, and so on. 
The obtained hierarchy is nothing but the Balitsky hierarchy [18] (see also 
Ref. [46,47]) restricted to dipoles. We refer the reader to [48,49] (see also 
Ref. [50]) for the detailed relationship of this equation to the B-JIMWLK 
formalism. Note in particular that the factorized correlators (10) is a solution 
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x o x o 



S(y,(Xo ,x 2 )) 



S(y,(^ ,x 3 )) 
S(y,(x 3 ,Xj )) 



X 2 



X, 



(a) 




Fig. 6. (a) Contribution to the B-JIMWLK equation for the 2-point correlator re- 
stricted to dipoles (x'2 is taken equal to X2 in this figure), (b) A graph that would 
also contribute to the 2-point correlator and that is missing in the B-JIMWLK 
formalism. 

of the whole hierarchy, and is actually a good approximation to the solution of 
the full B-JIMWLK equations. This statement was first made after the results 
of the numerical solution to the JIMWLK equation worked out in Ref. [51]. 



We may wonder why there are no terms involving one-point functions in the 
right handside of the previous equation. Actually, such terms would correspond 
to graphs like the one of Fig. 6b, in which, for example, two dipoles merge. 
They are expected to occur if saturation is properly taken into account. While 
the restriction of the Balitsky equation to dipoles does not drastically change 
the solution for the scattering amplitudes, such terms would instead have a 
large effect, as we shall discover in the following. In the next section, we will 
explain why such terms are actually required for physical reasons. 



2.1.4 Saturation 

The BK equation may be well-suited for the ideal case in which the target 
is a nucleus made of an infinity of independent nucleons. But it is not quite 
relevant to describe the scattering of more elementary objects such as two 
dipoles (or two virtual photons, to be more physical). 

Indeed, following Chen and Mueller [52] (see also Ref. [53]), let us consider 
dipole-dipole scattering in the center-of-mass frame, where the rapidity evo- 
lution is equally shared between the projectile and the target (see Fig. 7a). 
Then at the time of the interaction, the targets are dipoles that stem from the 
branching of a unique primordial dipole. Obviously, the assumption of statis- 
tical independence of the targets, which was needed for the factorization (10) 
to hold, is no longer justified. 

So far, we have seen that nonlinear effects which go beyond the factorization 
formula (7) are necessary to preserve unitarity as soon as n ~ l/a 2 s . This 
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came out of an analysis of Eq. (7) in the restframe of the target. The rapidity 
2/bfkl at which the system reaches this number of dipoles and hence at which 
the BFKL approach breaks down may be found from the form of the typical 
growth of n with y, namely n(y) ~ e ay . Parametrically, 

2/bfkl ~ -log^-. (13) 

a aj 

Now we may go to the center-of-mass frame, where Eq. (7) with y' = y/2 would 
describe the amplitude in the absence of nonlinear effects. There, the typical 
number of dipoles in the projectile and in the target are well below 1/a 2 : 
^(|/bfkl/2) ~ l/a s . We actually see that the evolution of the dipoles in each 
of theses systems remains linear until y = 2?/bfkl- In that rapidity interval, 
nonlinear effects consist in the simultaneous scatterings of several dipoles from 
the target and the projectile but the evolution of n still obeys the BFKL 
equation (see Fig. 7a). Now, performing a boost to the projectile restframe, the 
evolution goes into the target. Formula (3) should then apply for the amplitude 
A. But if the evolution of the target were kept linear, then the amplitude 
would not be unitarity because the number of dipoles would be larger than 
1/ct 2 . Hence, through some nonlinear mechanism, which was represented by 
multiple scatterings between linearly evolving objets in the center-of-mass 
frame, the dipole number density has to be kept effectively lower than 1/a 2 
in order to preserve unitarity (see Fig. 7b). This is called parton saturation. 
The precise saturation mechanism has not been formulated in QCD. It could 
be dipole recombinations due to gluon fusion, multiple scatterings inside the 
target which slow down the production of new dipoles [54] (as in Fig. 7b), 
"dipole swing" as was proposed more recently [55,56], or any other mechanism. 

Hence, unitarity of the scattering amplitudes together with boost-invariance 
seem to require the saturation of the density of partons. 

A pedagogical review of saturation and the discussion of the relationship be- 
tween saturation and unitarity may be found in Ref. [57]. Original papers 
include Refs. [58,59]. 

2.1.5 The Pomeron language 

So far, we have presented in detail a s-channel picture of hadronic interac- 
tions, and it is in this formalism that we will understand most easily the 
link with reaction-diffusion processes. In the s-channel formulation, all the 
QCD evolution happens in the form of quantum fluctuations of the inter- 
acting hadrons. However, a picture maybe more familiar to the reader is a 
t-channel picture, where the rapidity evolution is put in the t-channel, while 
the projectile and target are in their bare states. This picture directly stems 
from the usual Lorentz-invariant formulation of quantum field theory, while 
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(a) 



(b) 



Fig. 7. (a) Scattering in the dipole model in the center-of-mass frame. The evolution 
is shared between the target and the probe. The amplitude is unitarized through 
the multiple scatterings occuring between the two evolved wave functions, (b) Boost 
of the previous graph to the restframe of the projectile. There is now twice as much 
evolution in the target and the nonlinear effects should occur inside its wavefunction, 
in the course of the evolution. They may take the form of "internal" rescatterings 
(as depicted), or dipole merging... 





diagrams 

Fig. 8. The BFKL Pomeron is a sum of i-channel gluon Feynman diagrams. 

the dipole picture (or the parton model) is derived in the framework of time- 
ordered perturbation theory. 



Classes of Feynam diagrams can be grouped into "Pomerons" (or Reggeized 
gluons, see Fig. 8), in terms of which scattering processes may be analyzed. (A 
pedagogical review on how to derive the BFKL equation in such a formalism 
is available from Ref. [60]). 

An effective action containing Pomeron fields and vertices may be constructed. 
In these terms, the s-channel diagrams of Figs. 5,7a may be translated in 
terms of the diagrams of Fig. 9. The effective action formalism was initially 
developped in Refs. [61,62,63]. More recently, there has been some progress 
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(a) 



(b) 



Fig. 9. (a) Example of a diagram contributing to the BK equation in the i-channel 
representation (see Fig. 5). The dashed lines represent Pomerons. The rapidity is 
proportional to the length of the Pomeron lines in the t-channel. (b) Pomeron rep- 
resentation of a class of diagrams to which Fig. 7a belongs. 

in the definition of the effective action [64], some of it with the help of the 
correspondence with statistical physics processes [65,66]. 

We will not expand on this formulation in the present review, because it 
is difficult to see the analogy with statistical processes in this framework. 
A s-channel picture is much more natural. However, a full solution of high- 
energy QCD may require to go back to that kind of calculation and compute 
accurately the 1 — > n Pomeron vertices. This program was formulated some 
time ago [67,68], and there is continuing progress in this direction (see e.g. 



2.2 Analogy with reaction- diffusion processes 

We are now in position to draw the relationship between high-energy QCD and 
reaction-diffusion processes. In the first section below, we will show that the 
BK equation is, in some limit, an equation that also appears in the context 
of statistical physics. Second, we will exhibit a particular reaction-diffusion 
model, and show in the final section how this model is related in a more 
general way to scattering in QCD. 

2.2.1 The BK equation and the FKPP equation 

Let us first show at the technical level that under some well-controlled approx- 
imations, the BK equation (11) may be mapped exactly to a parabolic nonlin- 



[69,70]). 
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ear partial differential equation. This observation was first made in Ref. [38]. 

To simplify, we will look for impact-parameter independent solutions: A(y, (x , Xi)) 
is supposed to depend on y and x 01 only, not on x + x±. We switch to mo- 
mentum space through the Fourier transformation 

A(y,k) = J ^/ kxoi My,xoi). (14) 

This transformation greatly simplifies the BK equation [24,25]. It now reads 

da y A(y, k) = X (-d, ogk 2)A(y, k) - A 2 (y, k). (15) 

The first term in the right handside, which is a linear term, is actually an 
integral kernel, obtained by Fourier transformation of the BFKL kernel (first 
three terms in the right handside of Eq. (11)). It is most easily expressed 
in Mellin space: k~ 21 is the set of its eigenfunctions, with the corresponding 
eigenvalues 

x( 7 ) = 2^(1) -#y)- ^(1-7)- (16) 
This kernel may be expanded around some real 7 = 70, fixed between and 
1. Keeping the terms up to 0((j — 70) 2 ) is the well-known diffusive approx- 
imation, which is a good approximation for large rapidities. Introducing the 
notations Xo = x(7o), Xo = x'(To) and x'o = x"(7o), the BK equation reads, 
within this approximation 

d &y A = 4dl gk *A + (70X0 - Xo)dio Sk >A + ( X o - I0X0 + 2 ¥ 1 )^ - A 2 . (17) 
Through some linear change of variable (ay, log A; 2 ) — > (t,x), 

t 



ay = - 

Xo — I0X0 

log/c 2 = 



Xo ToXo ~ Xo , 



(18) 



\ 2(xo - I0X0) + 7oXo' Xo ~ I0X0 + 



one may get rid of the first-order partial derivative in the right handside. We 
then find that the new function 

Xo ~ 7oXo + ^ 



obeys the equation 

du(t,x) d 2 u(t,x) 
dt dx 2 



+ u(t,x) -u 2 (t,x), (20) 



which is the Fisher [36] and Kolmogorov-Petrovsky-Piscounov [37] (FKPP) 
equation. This equation was first written down as a model for gene propagation 
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in a population in the large population size limit. But it turns out to apply 
directly or indirectly to many different physical situations, such as reaction- 
diffusion processes, but also directed percolation, and even mean-field spin 
glasses [71]. A recent comprehensive review on the known mathematics and 
the phenomenological implications of the FKPP equation can be found in 
Ref. [72]. 

As a side remark, we note that if 70 is chosen such as xilo) = lox'ilo), then the 
mapping drastically simplifies. Actually, this choice has a physical meaning, 
as we will discover in Sec. 4 when we try and solve the BK equation. 

Beyond the exact mapping (20) between an approximate form of the BK and 
the FKPP equations, the full BK equation is said to be in the universality 
class of the FKPP equation. All equations in this universality class share 
some common properties, as will be understood below. The exact form of 
the equation is unessential. As a matter of fact, recently, it has been checked 
explicitely that the BFKL equation with next-to-leading order contributions 
to the linear evolution kernel (but keeping the QCD coupling fixed) is also 
in the same universality class. A mapping to a partial differential equation 
(which involves higher-order derivatives in the rapidity variable) was exhibited 
[73] . What defines physically the universality class of the FKPP equation is 
a branching diffusion process with some saturation mechanism. The details 
seem unimportant. 

In the next section, we shall give a concrete example of a reaction-diffusion 
process: We will see how the FKPP equation appears as a fluctuationless (or 
"mean-field") limit of some stochastic reaction-diffusion process. In Ref. [38], 
it had not been realized that the analogy of QCD with such processes is in 
fact much deeper than the formal mapping between the BK equation and the 
FKPP equation that we have just outlined. But this is actually the 
we shall shortly argue. 

2.2.2 Reaction- diffusion processes: An example 

We consider a reaction-diffusion model, which was introduced in Ref. [74]. 
Particles are evolving in discrete time on a one-dimensional lattice. At each 
timestep, a particle may jump to the nearest position on the left or on the 
right with respective probabilities pi and p r , and may split into two particles 
with probability A. We also allow that each of the n(t, x) particles on site x 
at time t to die with probability Xn(t,x)/N. 

From these rules, we may guess what a realization of this evolution may look 
like at large times. The particles branch and diffuse (they undergo a linear 
evolution) until their number n becomes of the order of N, at which point the 
probability that they "die" starts to be sizable, in such a way that their number 
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n(t,x) n(t+dt,x) 




Fig. 10. Picture of a realization of the system of particles at two successive times. 
In the bins in which the number of particles is of order N, some particles disappear, 
others are created by splittings, but overall the number of particles is conserved 
up to fluctuations of order y/N. In the bins in which n is small compared to N, 
the dynamics is driven by branching diffusion. As a result, n(t, x) looks like a noisy 
wave front moving to the right. 

never exceeds N by a large amount, on any site. But if the initial condition 
is spread on a finite number of lattice sites, the linear branching- diffusion 
process may always proceed towards larger values of \x\, where there were no 
particles in the beginning of the evolution. Hence a realization will look like 
a front connecting an ensemble of lattice sites where a quasi-stationary state 
in which the number of particles is N (up to fluctuations) has been reached, 
to an ensemble of empty sites (towards \x\ — > oo). This front moves with 
time as the branching diffusion process proceeds. The position of the front 
X(t) may be denned in different ways, leading asymptotically to equivalent 
determinations, up to a constant. For example, one may define X(t) as the 
rightmost bin in which there are more than N/2 particles, or, alternatively, 
as the total number of particles in the realization whose positions are greater 
than 0, scaled by 1/N. A realization and its time evolution is sketched in 
Fig. 10. 

Between times t and t + At, ni(t,x) particles out of n(t,x) move to the left 
and n r (t,x) of them move to the right. Furthermore, n + (t,x) particles are 
replaced by their two offspring at x, and n-(t,x) particles disappear. Hence 
the total variation in the number of particles on site x reads 

n(t + At, x) — n(t, x) = —ni(t, x) — n r (t, x) — n_(t, x) 

+ n + (t, x) + ni(t, x + Ax) + n r (t, x — Ax). (21a) 

The numbers describing a timestep at position x have a multinomial distribu- 
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tion: 



n< 

P({n h n r ,n + ,n-}) = — — : — - — -— -p? l p nr 

n|!n r !n + !n_!An! 



r 

An 



x X n +(Xn/N) n -(l-pi-p r -X-Xn/N) An , (21b) 

where An = n — ni — n r — n + — n_ , and all quantities in the previous equation 
are understood at site x and time t. The evolution of u = n/N is obviously 
stochastic. One could write the following equation: 

u(t + At, x) = (u(t+At, x)) + ^f(u 2 (t+At,x)) - {u(t+At,x)) 2 u(t + At, x) 

(22) 

where the averages are understood over the time step that takes the system 
from t to t + At. They are conditioned to the value of u at time t. v is a noise, 
i.e. a random function. The equation was written in such a way that it has 
zero mean and unit variance. Note that the noise is updated at time t + At in 
this equation. 

One can compute the mean evolution of u = n/N in one step of time which 
appears in the right handside of Eq. (22) from Eq. (21). It reads 

(u(t+At, x)\{u(t, x)})=u(t, x)+pi[u(t, x+Ax) —u{t, x)] 

+p r [u(t, x-Ax)-u(t, x)]+Xu(t, x) [l-u(t, x)]. (23) 

The mean evolution of the variance of u that appears in Eq. (22) may also 
be computed. The precise form of the result is more complicated, but roughly 
speaking, the variance of u after evolution over a unit of time is of the order 
of u/N for small u ~ 1/N. This is related to the fact that the noise has a 
statistical origin: Having n particles on the average in a system means that 
each realization typically consists in n ± y/n particles. 

When N is infinitely large, one can replace the w's in Eq. (23) by their averages: 
This would be a mean field approximation. Obviously, the noise term drops 
out, and the equation becomes deterministic. Note that if we appropriately 
take the limits Ax — > and At — > 0, setting 

At 

A = Af, PR = PL = j£^, (24) 

the obtained mean-field equation is nothing but the FKPP equation (20). For 
the numerical simulations of this model that we will perform in Sec. 4, we 
will keep At and Ax finite, which is usually more convenient for computer 
implementation. 

Thus we have seen that the evolution of reaction-diffusion systems is governed 
by a stochastic equation (22) whose continuous limit (At — > 0, Ax — > 0) 
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and mean-field limit (N ^> 1) is a partial differential equation of the form of 
(exactly actually, in our simple case study) the FKPP equation. We shall now 
argue that partons in high-energy QCD form such a system. 



2.2.3 Universality class of high- energy QCD 

Let us come back to the QCD dipole model. We have seen that rapidity 
evolution of the hadron wavefunctions proceeds through a branching diffusion 
process of dipoles. Let us denote by T(y, r) the scattering amplitude of the 
probe dipole off one particular realization of the target at rapidity y and at a 
given impact parameter. This means that we imagine for a while that we may 
freeze the target in one particular realization after the rapidity evolution y, 
and probe the latter with projectiles of all possible sizes. Of course, this is not 
doable in an actual experiment, not even in principle. But it is very important 
for the statistical picture to go through such a "gedanken observable". The 
amplitude A, which is related to the measurable total cross section, is nothing 
but the average of T over all possible realizations of the fluctuations of the 
target, namely 

A(y,r) = (T(y,r)). (25) 

The branching diffusion of the dipoles essentially occurs in the log(l/r 2 ) vari- 
able. The scattering amplitude is roughly equal to the number of dipoles in a 
given bin of (logarithmic) dipole size, multiplied by a 2 . From unitarity argu- 
ments and consistency with boost-invariance, we have seen that the branching 
diffusion process should (at least) slow down in a given bin as soon as the num- 
ber of objects in that very bin is of the order of N = l/a 2 , in such a way that 
effectively, the number of dipoles in each bin is limited to N. A typical re- 
alization of T is sketched in Fig. 11. As in the case of the reaction-diffusion 
process, from similar arguments, it necessarily looks like a front. The position 
of the front, defined to be the value r s of r for which T is equal to some fixed 
number, say |, is related to the saturation scale defined in the Introduction: 

r s = i/Q s {y)- 

We now see that there is a very close analogy between what we are describing 
for QCD here and the model that we were introducing in the previous section. 
So in particular, one might be able to formulate interaction processes in QCD 
with the help of a stochastic nonlinear evolution equation for the "gedanken" 
amplitude T. We already know the mean-field limit that one should get, when 
N is very large: This is the BK equation, as was rigorously proven above. Thus 
we know the equivalent of the term (u(t + At, x)) in Eq. (22). The noise term 
is not known, but since it is of statistical origin, we know that it must be of 
the order of the square root of the number of dipoles normalized to N, that 
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log(r/r 2 ) 



Fig. 11. Sketch of the scattering amplitude T of a dipole of size r off a frozen 
partonic configuration. The small lines on the axis denote the dipoles ordered by 
their logarithmic sizes. Up to fluctuations, T looks like a wave front. 



is to say, of order JT/N. We may write an equation of the form 



d ay T(y, k) = x(-di ogk 2)T(y, k) - T 2 (y, k) + a s ^2T(y, k) v{y, k), (26) 

where v is a noise, uncorrelated in rapidity and transverse momentum, with 
zero mean and unit variance. (The factor of 2 under the square root is essen- 
tially arbitrary). This equation is to be compared to the following one: 



d t u(t,x) = d 2 x u{t,x) +u(t,x) -u 2 (t,x) + ^/ 2m (^ x ) u (t,x), (27) 

which is the so-called "Reggeon field theory" equation when the noise v is 
exactly a normal Gaussian white noise, that is to say, of zero mean and whose 
non- vanishing cumulant reads 

(u(t, x)u(t', x')) = 5(t - t')5(x - x'). (28) 

It is a stochastic extension of Eq. (20). If the noise term were of the form 



I2u(t,x)(l-u(t,x)) 



N 



u(t, x) 



(29) 



instead, then this equation would be what is usually referred to as the stochas- 
tic Fisher-Kolmogorov-Petrovsky-Piscounov equation. The sFKPP equation 
and the physics that it represents is reviewed in Ref. [75]. 

Taking averages over events converts this equation into a hierarchy of cou- 
pled equations, which has a lot in common in its structure with the Balitsky 
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Reaction-diffusion 



QCD 



Occupation fraction u(t, x) 



Scattering amplitude for the probe off a 
frozen realization of the target T(k, y) 

Physical scattering amplitude A = (T) 



Average occupation fraction (u(t,x)} 



Space variable x 



\og{k 2 /K 2 ) or log(l/r 2 A 2 ) 



Time variable t 




Average maximum density of particles N 



Position of the front X(t) 

Branching-diffusion kernel u{— d x ) 
(uj(-d x ) = dl + 1 in the FKPP case) 



BFKL kernel xi-diogtf) 

or its equivalent in coordinate space 



Table 1 



Dictionary between QCD and the reaction-diffusion model for the main physical 
quantities. A is a typical hadronic scale. 

hierarchy (12). (Actually, there are some extra terms compared to the Balit- 
sky hierarchy, which were first found from the analogy with reaction-diffusion 
processes, and which precisely represent nonlinear effects inside the wavefunc- 
tions. A detailed study may be found in Ref. [76]). We will perform explicit 
calculations in this spirit within simpler models in Sec. 3 below. 

Based on these considerations, we may establish a dictionary between QCD 
and reaction-diffusion processes. The correspondence is summarized in Tab. 1. 

The mechanism for saturation of the parton densities (i.e. of the dipole number 
density) is not known for sure in QCD. There are also important differences 
between the reaction-diffusion model introduced above and QCD that lie in the 
"counting rule" of the particles (provided by the form of T el in the QCD case, 
see Eq. 2). But from the general analysis of processes described by equations 
in the universality class of the stochastic FKPP equation and the underlying 
evolution mechanisms presented in Sec. 4, we will understand that most of the 
observables have universal properties in appropriate limits, which do not de- 
pend on the details of the mechanism at work. We draw the reader's attention 
to Refs. [77,78], where a precise stochastic equation was searched for in QCD. 
Some of the problems one may face with the use and the very interpretation 
of such equations were studied in Ref. [79]. 

The way in which we view high energy QCD is actually not particularly orig- 
inal: It is nothing but the QCD dipole model, which was implemented numer- 
ically in the form of a Monte Carlo event generator by Salam [80,81,82] (see 
also [55] for another more recent implementation). He also devised and im- 
plemented a saturation mechanism [54] that went beyond the original dipole 
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model pictured in Fig. 7a, but which is necessary, as was argued before. 



Before discussing more deeply the physical content of equations of the form 
of Eq. (26), we shall first study a model in which spatial dimensions are left 
out, that we will be able to formulate in different ways. 



3 Zero-dimensional model 

In the previous section, we have understood that scattering at high energy 
in QCD may be viewed as a branching-diffusion process supplemented by 
a saturation mechanism. We have exhibited a simple toy model with these 
characteristics, whose dynamics is represented by an equation of the type (27). 

Unfortunately, even that toy model is too difficult to solve analytically. We 
shall study a still simplified model, where there is no diffusion mechanism: Re- 
alizations are completely specified by the number of particles that the system 
contains at a given time. Of course, in this case, a saturation scale cannot be 
defined, which limits the relevance of this model for QCD. However, we will 
be able to formulate this model in many different ways, and to draw parallels 
with QCD. 

We start by defining precisely the model. Then, two approaches to the com- 
putation of the moments of the number of particles are presented. The first 
set of methods relies on field theory (Sec. 3.2). The second method relies on a 
statistical approach (Sec. 3.3) and will be extended in a phenomenological way 
to models with a spatial dimension in Sec. 4. We shall then draw the relation 
to a scattering-like formulation (Sec. 3.4). Finally (Sec. 3.5), some variants of 
the basic model are reviewed. 

3. 1 Definition 

Let us consider a simple model in which the system is characterized by its 
number n t of particles at each time t. Between times t and t + dt, each particle 
has a probability dt to split in two particles. For each pair of particles, there is 
a probability dt/N that they merge into one. We may summarize these rules 
in the following form: 



n t + l proba n t dt 



n t — l proba 



n t (n t -l)dt 



nt+dt 



N 



(30) 



n t proba l—n t dt — 



n t {n t -l)dt 
N 
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From this, one can easily derive an equation for the time evolution of the 
probability P(n, t) of having exactly n particles in the system at time t: 



^(n, t) = (n- l)P(n -l,t) + ^^P(n + l,t)-(n + ^j^j P(n, t). 

(31) 

This is the master equation for the Markovian process under consideration. 
The two first terms with a positive sign represent the process of going from 
one state containing n particles to an adjacent one containing n + 1 or n — 1 
particles respectively, while the last term simply corrects the probability to 
keep it unitary. 

By multiplying both sides of this equation by n and summing over n, we get 
an evolution equation for the average number of particles (n t ): 

^ = <n t >-l<n t (n t -l)> (32) 

Unfortunately, this equation is not closed, and one would have to establish an 
equation for (n t (n t — 1)), which would involve 3-point correlators of n t , and 
so on, ending up with an infinite hierarchy of equations, exactly like in Sec. 2 
for QCD (see Eq. (12)). 

This illustrates the difficulties one has to face before one can get an analytical 
expression for (nt), even in such a simple model. 



3.2 "Field theory" approach 



In the next subsections, we will follow different routes to get analytical results 
on the moments of the number of particles in the system at a given time t. The 
first one will be similar to the s-channel picture of QCD (see Sec. 2), since it 
will consist in computing the time (rapidity in QCD) evolution of realizations 
of the system. The second one will be closer to the t-channel picture of QCD. 
We will see how "Pomerons" may appear in these simple systems. We will then 
examine a formulation in terms of a stochastic nonlinear partial differential 
equation, which is nothing but the sFKPP equation in which the space variable 
(x) has been discarded. 



3.2.1 Particle Fock states and their weights 

Statistical problems were first formulated as field theories by Doi [83] and Peliti 
[84]. Different authors have used these methods (see Ref. [85] for a review). 
We shall start by following the presentation given in Ref. [86]. 
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We would like to interpret the master equation (31) as a quasi-Hamiltonian 
evolution equation of the type of the ones that appear in quantum mechanics. 
To this aim, we need to introduce the basis of states \n) of fixed number n 
of particles. We define the ladder operators a and a* by their action on these 
states: 

a\n) = n\n — 1), a)\n) = \n+ \) (33) 
and which obey the commutation relation 

[a,a f ] = l. (34) 

The n-particle state may be constructed from the vacuum (zero-particle) state 
by repeated application of the ladder operator: 



\n) = (a*y\0). (35) 



The normalization is not standard with respect to what is usually taken in 
quantum mechanics. In particular, the orthogonal basis \n) is normalized in 
such a way that (m\n) = n\5 m , n . This implies that the completeness relation 
reads 

E^|n)(n| = l. (36) 

We also introduce the state vector of the system at a time t as a sum over all 
possible Fock states weighted by their probabilities: 

\m=Y, P (n,t)\n). (37) 

n 

It is straightforward to see that the master equation (31) is then mapped to 
the Schrodinger-type equation 

J^(f)> = -H\<f>(t)), (38) 
where Tt is the "Hamiltonian" operator 

H = (l-o t )o t o--(l-o t )aV. (39) 

The first term represents the splitting of particles, while the second one, pro- 
portional to 1/N, represents the recombination. We may rewrite Ti as 

H = H + Hi, (40) 

where 

Ho = a f a (41) 

is the "free" Hamiltonian whose eigenstates are the Fock states. We now go 
to the interaction picture by introducing the time-dependent Hamiltonian 

Hi(t) = e Hot Hie- not (42) 
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and the states \<f>)i = e Hot \4>). The solution of the evolution reads 

|0) / = Texp (- fdt'Hiit'))^)! 

t t («) 

= \<h)i-l dt"H I (t > )\ ( f>o)i+ [ dt' [ dfHi(i!)Hi{f)\<h)i + '-- 

JO Jo JO 

We may then compute the weights of the successive Fock states by applying 
this formula. Let us show how it works in detail by computing the state of a 
single particle evolved from time to time t, in the limit N = oo in which 
there are no recombinations. We follow the usual method to deal with such 
problems in field theory. We insert repeatedly complete basis of eigenstates of 
Ho into Eq. (43), namely 

10)/ = |1> " fdt'Y. -^rlnxXmlW^OIl) + • ■ • (44) 

JO ni Tl\. 



(We have kept the first two terms in Eq. (43) explicitely) . Using the expression 
for Hi(t) as a function of Ho and Hi, together with the knowledge that the 
Fock states are eigenstates of H , we get 

10) = e"*|l> - E^" 1 * f dt'e^ t '- t ' — \n 1 )(n 1 \H 1 \l) + ■■■ (45) 



Inserting the expression for Hi, one sees that in the infinite- N limit, there 
is only one possible transition, 1 — > 2. Performing the integration over t' and 
computing in the same manner the higher orders, one finally gets the expansion 

|0) = e -*|l) + e -*(l - e -*)|2) + • • • + e~\\ - e-*)"-» + • • • (46) 

from which one can read the probabilities of the successive Fock states. This 
expansion is similar to the expansion in dipole Fock states introduced in Sec. 2: 
The n-particle states correspond to n-dipole states in QCD, and their weights 
are computed by applying successive splittings to the system, whose rates are 
given by Eq. (1). (They are just unity in the case of the zero-dimensional 
model.) 

We see that this method is well-suited to compute the probabilities of the 
lowest-lying Fock-states, and their successive corrections at finite N . But in 
general we are rather interested in averages such as (n k ), for which the weights 
of all Fock states are needed. We will develop a slightly different (but equiv- 
alent) formalism below, that will enable us to get these averages in a much 
more straightforward way. 
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3.2.2 Pomeron field theory 

Let us introduce the generating function of the factorial moments of the dis- 
tribution of the number of particles 

Z(z,t)=J2(l + z) n P(n,t). (47) 

n 

The evolution equation obeyed by Z can easily be derived from the master 
equation (31): 

dz m \ ( dz 1 d2z \ 

m= z{1 + z) [^-N^)- (48) 

We may represent this equation in a second-quantized formalism by introduc- 
ing the operators 

6 f = z, b=^-=z (49) 
oz 

acting on the set of states \Z) consisting in the analytic functions of z. Then 
we may write 

f = -*Z, (50) 

where 

H p = Hi + Hi with Hi = -b% H\ = -bWb + lftt(l + tf)b\ (51) 
A basis for the states is 

\k)=z\ (k\ = z k (52) 
which is orthogonal with respect to the scalar product 

(ZAZ 2 ) = J ^ e ^ 2 Z 1 (z,z)Z 2 (z,z), (53) 

and obeys the normalization condition (k\l) = k\5k,i- We shall call these states 
"/c-Pomeron" states, by analogy with high-energy QCD. We may apply exactly 
the same formalism as before, since the operators b, have the same properties 
as the a, a) . 

From the definition of the scalar product, it is not difficult to see that the k- 
th factorial moment of n may be obtained by a mere contraction of the state 
vector \Z), computed by solving the Hamiltonian evolution, with a A;-Pomeron 
state. The following identity holds: 



<*!*> = (54) 



where the average in the right handside goes over the realizations of the system. 
As for the initial condition, starting the evolution with one particle means 
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t 

(a) (b) (c) (d) 

Fig. 12. Propagator and vertices for the Pomeron field theory. Time flows from the 
top to the bottom. 

taking as an initial condition the superposition |0) + |1) of zero- and one- 
Pomeron states respectively. The zero-Pomeron state does not contribute to 
the evolution, hence a one-Pomeron state is like a one-particle state. 

In order to simplify the systematic computation of these moments, we may use 
a diagrammatic method and establish Feynman rules. To this aim, we write 
the contribution of the graphs with /-vertices (corresponding to the term of 
order / in the expansion of Eq. (43)), starting with a one-Pomeron state: 

(k\z) d (-i) L fdh f 1 dt 2 --- f" 1 dt t £ (k^—MTfM^) ■ ■ • J-(m|w?|i). 

JO Jo Jo ni7^,n« n l ] n l ! 

(55) 

Each matrix element that appears in this equation is associated to a vertex, 
and propagators connect these vertices. We read off the expression for the 
Hamiltonian (51) that there is one propagator and three vertices in the theory: 
one splitting (1 — > 2), one recombination (2 — > 1) and a 2 — > 2 elastic diffusion 
vertices. 



The method to compute the 1 to k Pomeron transition amplitude is standard. 
First, one draws all possible diagrams for this transition that contain / ver- 
tices, including all possible permutations. (Note that a splitting may occur 
in k different ways, if k is the number of Pomerons before the splitting; A 
recombination instead may occur in k(k — l)/2 ways). Then, the propagators 
(Fig. 12a) are replaced by 

(l|e- m o|l) = e *, (56) 

(where t is the time interval that they span) in such a way that the n-Pomeron 
state propagates as (n\e~ tn o\n) = e nt . Intermediate times are integrated out. 
As for the vertices (Figs. 12b-d), the following factors have to be applied: 

2 2 
l-»2: -1; 2->l: — ; 2^2: — . (57) 

In addition, there is a (_i)# vcrtlces factor. Finally, an overall k\ factor leads to 
the expression for the factorial moment (n t (n t — 1) • • • (n t — k + 1)). 
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(k! such diagrams) 

Fig. 13. Diagrams contributing to the one Pomeron — > /c-Pomeron transition, which 
gives the moments (nt(nt — 1) ■ ■ ■ (n< — k + 1)) at leading order in a 1/iV expansion. 

The lowest-order diagram for the average particle number, consisting in a sim- 
ple propagator, reads (n t ) = e*. We now understand that this method leads to 
a more straightfoward computation of the moments of the number of particles 
than the one based on the computation of the probabilities of successive Fock 
states, for a single Pomeron already resums an infinity of particle Fock states. 
The Pomeron in this case is exactly like the BFKL Pomeron introduced in 
Sec. 2, which leads to an exponential increase of the scattering amplitudes 
with the rapidity (Eq. (56)). 

We now move on to the computation of higher-order diagrams. First, let us 
recover simple results by taking the infinite- N limit. We consider the diagrams 
in Fig. 13, which are the only ones that survive for N = oo in the evaluation 
of the moment (n t (n t — 1) • • • (n t — k + 1)). Using the Feynman rules, we get 
for each individual diagram 

(-l) fc x (-l) fe x e kt f'dhe-' 1 C dt 2 e~ t2 ■ ■ ■ f dt k e~ tk = — ( 1 - e-*)^ 1 . 

JO Jt x Jtk-i 

(58) 

There are k\ such diagrams (corresponding to all possible permutations of the 
Pomerons), and there is an extra overall k\ factor to be added in order to get 
the relevant factorial moment: 

(n t {n t - 1) • • • (nt - k + 1)) = k\e kt (l - . (59) 



Next, we would like to perform the computation of the one- Pomeron — > one- 
Pomeron transition (which provides the value of (n t )) within the full theory, 
including the recombinations. Some of the lowest-order diagrams are shown 
in Fig. 14. A straightforward application of the Feynman rules edicted above 
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(a) 



(b) 



(c) 



(d) 



(e) 



Fig. 14. Diagrams up to order l/N 2 contributing to the average of the number of 
particles in the system after an evolution over the time interval t. 

leads to the following results for the graphs that are depicted in Fig. 14: 



( n *) tree, Fig. 14a 




\ n t)\l loop, Fig. 14b 


2t 


= -2! — 

N 


( n t) 2 loops, Fig. 14c 




= 3!— ( 

N 2 V 






\ n t)\2 loops, Fig. 14d 




( n t) 2 loops, Fig. 14c 







(1- 



1 + i 



3!— (l + 4e-*(l-t)- e - 2 '(2t + 5)) 



(60) 



We may classify these different contributions according to their order in e t /N: 
The leading terms for large t and e t /N ~ 1 are always of the form A^(e'/^) 1+ * loops - 
It turns out that we may compute easily these dominant terms at any number 
of loops. They stem from the graphs in which all splittings occur before all 
recombinations. These terms build up a series that reads 



(nt) = Et- 1 )^' 



k=l 



(61) 



This series is factorially divergent, but is easy to resum with the help of the 
Borel transformation. Indeed, using the identity 



k\ 



! = / dbb k e-\ 
Jo 



(62) 



then exchanging the integration over b and the sum over the number of 
Pomerons k, one gets 



r+oo 

(nt) = iVV* / db 
Jo 



-Ne-'b 



= N(l-Ne Ne ~ t T{0,Ne- t )), (63) 
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where T is the incomplete Gamma function. 

This result was obtained for the first time using a diagrammatic method in 
Ref. [87]. The authors of that paper also computed the next-to- leading order, 
that is to say, the terms of relative order 1/N after the resummation has been 
performed. The equivalent of the diffractive processes known in QCD were 
also investigated by the same authors in Ref. [88]. More results were obtained 
on that kind of models by another group in Ref. [89,90], using different tech- 
niques, which go beyond the perturbative approach. Remarkably, the latter 
calculations can be applied to some extent to QCD [91,92]. 



3.2.3 Stochastic evolution equations 

The model may also be formulated in the form of a stochastic evolution equa- 
tion for the number of particles nt it contains at each time t. The most straight- 
forward way of doing this would be to first compute the mean and variance of 
n t+ dt given n t , with the help of the master equation (31). This would enable 
one to write the time evolution of n t in terms of a drift and of a noise of zero 
mean and normalized variance, namely: 



dn t n t (n t -l) / n t (n t - 1) 

where v is such that [v t ) = and (ytVf) = 5(t — t'). This equation is similar to 
Eqs. (26) and (27), except for it does not have a spatial dimension where some 
diffusion could take place. The noise term is of order s/n, as it should according 
to the argumentation of Sec. 2. Note that the distribution of v depends on n t 
and is not a Gaussian. This last point is easy to understand: The evolution 
of v t is intrinsically discontinuous, since it stems from a rescaling of n t , which 
is an integer at all times. A Brownian evolution (i.e. with a Gaussian noise) 
would necessarily be continuous. For completeness, we write the statistics of 
vt+dt, which is easy to derive from the evolution of n: 



vt+dt = { -£ proba l-n t dt-^p^dt (65) 



_L 


A 


a dt 






A 




(7 


1 


A 


a dt 




and a 





proba n t dt 

N 



proba ^^dt, 



where A = n t — nt ^ ^ and a = yn t + nt ^ There are jumps induced by 
the terms proportional to 1/dt. 

This formulation is not of great interest, neither for analytical calculations nor 
for numerical simulations, since it is much easier to just implement the rules 
that define the model in the first place (Eq. (30)) in the form of a Monte Carlo 
event generator. 



35 



There is a better way to arrive at a stochastic evolution equation for this 
model, although it is a bit more abstract. (It is actually equivalent to the 
Pomeron field theory formulated before.) Instead of following states with a 
definite number of particles like above, we may introduce coherent states 



\z) 



-z+za' 



|0>, 



(66) 



where z is a complex number. For real positive values of z, the state \z) is 
nothing but a Poissonian state, which is a superposition of |/c)-particle states, 
where the weight of each state follows the Poisson law of parameter z. For 
the simplicity of the argument, let us restrict ourselves to these states. By 
applying the Hamiltonian H (defined in Eq. (39)) to a Poissonian state \z t ), 
one gets a new state \4>t+dt)'- 



\<f>t+dt) = \zt) - dtH\z t ) 



(67) 



Of course, that new state is not itself at Poissonian state in general, but may 
be written as a superposition of such states. One writes 



\M = Jdzf(z)\ z } = fdzWE 



e —n. 



(68) 



The idea is to interpret the weight function f(z) as the probability to observe 
a given Poissonian state \z). Hence the evolution is viewed as a stochastic path 



Zt-dt 



z t 



z t+dt — * Zt+2dt 



(69) 



with well-defined transition rates from one Poissonian state to the next one. 
Inserting the explicit expression for the Hamiltonian (39) and the decomposi- 
tion (68) in Eq. (67), one gets for each Fock state \n) 



1 



dze- z f{z)- 



e -«fL 
n! 



dt e~ 



zt 



Zt 



n-1 



n+1 



(n-1)! (n-2)! 




(n-2)! 



(70) 



Finally, this equation is easy to invert for f(z) by integrating over n with the 
weight / a l on g an appropriate contour in the complex plane. After 

some straightforward algebra, we get 



f{zt+dt) = S(zt+dt - z t ) +dt[zt 




zt 



%)5"(z t+dt -z t ) 



(71) 



This is a Gaussian centered at z t + dt(z t — %) of variance 2dt{z t -%). Intro 
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ducing a normal Gaussian noise v t which satisfies 



(v t ) = and (u t u t >) = 5(t-t'), 



(72) 



we may write 




(73) 



where the noise is taken at time t + dt, and hence, this equation is to be 
interpreted in the Ito sense. If z t =o is a real number between and N, then 
the equation keeps it in this range. But one may consider more general coherent 
states, with complex z t . 

This equation is suitable for numerical simulations: One may discretize the 
time in small steps At -C 1 in which case v t is distributed as 



p{y t ) 



l 



V2nAt 



cxp 



2At 



(74) 



(In many cases, one has to use more sophisticated methods, see e.g. Ref. [93]). 
Analytical manipulations of this equation using Ito's calculus are also quite 
easy. We are going to give an example of such a calculation below, avoiding 
unnecessary formalism. (We refer the reader to [93] for a textbook on a more 
mathematical treatment of stochastic processes.) 

We may transform the stochastic equation (73) to a hierarchy of equations for 
the factorial moments of the number of particles, using the relation 



(z k t ) = {n t {n t -l)---{n t -k+l))=nf\ 
First, let us write Eq. (73) in a discretized form: 



(75) 



Zt+dt — z t + dt [ z t 



z 



n) +dt 



\ 



y 2\ 



2U 



N 



I v t+ dt- 



(76) 



We then take the k-th power of the left and the right handside, and we average 
the result over realizations. Expanding in powers of dt for small dt, we get 



't+dt / 



(z k t ) + dtk(z k t - 



Jc+1 



TV 



+ dtk( z: 



k-l 



2U- 



N . 



(Vt+dt, 



+ dt 



2 



t+dt) 



+ 



(77) 



We have factorized the average over the noise over the time intervals [t, t + dt] 
and [0, t], since the noise v is uncorrelated in time. The term proportional to dt 
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vanishes thanks to the fact that vt+dt averages to zero. One may think that the 
next term could be neglected for it is apparently proportional to dt 2 . Actually, 
it gives a contribution of order dt, because for discretized t, (vt+dt) 

= 1/dt. 

The dots stand for terms of order dt 2 at least. Using Eq. (75) to identify the 
factorial moments of n, we eventually get 

^ = *(* w -^) + *<*- 1 >('^ I, -7r I (78) 

This equation is similar to the (modified) Balitsky hierarchy in high-energy 
QCD. Let us write explicitely the equations for the first two moments: 

d{n t (n t - 1)) 



dt 



2 (l - ^) (n t (n t - 1)) - ^(n t (n t - l)(n t - 2)) + 2(n t ). 

(79) 



We note the similarity in structure with Eq. (12), except for the term 2(n) in 
the right handside. This term stems precisely from the particle recombinations, 
and was absent in the B-JIMWLK formalism. 

Finally, let us mention that for a more rigorous and general derivation of this 
stochastic formulation, one may use a path integral formalism obtained from 
the Hamiltonian (39), see Ref. [85]. 



3.3 Statistical methods 



The field theory methods presented above provide a systematics to solve the 
evolution of the system to arbitrary orders in l/N, at least theoretically. (In 
practice, identifying and resumming the relevant diagrams becomes increas- 
ingly difficult). However, it would look quite unreasonable to get into such an 
involved formalism if one were only interested in computing the lowest order 
in a large- iV expansion. Indeed, as we shall demonstrate it below, in the case 
of this simple model, an intuitive and economical calculation leads to the right 
answer [94]. We work it out here because this line of reasoning is at the ba- 
sis of the solution to more complicated models, closer to QCD, that we shall 
address in the next section (Sec. 4). 

As before, we denote by n t the value of the number of particles in a given 
realization of the system. We further introduce pn{t) the distribution of the 
times at which the number of particles in the system reaches some given value 
n, and (n t \n, t) the conditional average number of particles at time t given that 
there were n particles in the system at time i. One may write the following 
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Fig. 15. [From Ref. [94]] Ten different realizations of the stochastic evolution of 
the zero-dimensional model (dotted lines; N = 5 x 10 3 ). All realizations look the 
same, up to a shift in time. They are all parallel to the solution to the mean-field 
equation (81) (dashed line). Note the significant difference between the latter and 
the average of the particle number over the realizations (full line). 

factorization formula: 



This formula holds exactly for any value of n. In particular, if iV is large 
enough, one may choose n such that 1 <C n <C N. 

Looking at a few realizations generated numerically (Fig. 15), one sees that 
the curves that represent n t look like the solution to the mean-field equation 
obtained by neglecting the noise term in Eq. (64), up to a translation of the 
origin of times by some random to. (The curves look also slightly noisy around 
the average trend, but the noise would still be much weaker for larger values 
of N.) This suggests that once there are enough particles in the system (for 
n t > n ^> 1), the evolution becomes essentially deterministic and in that 
stage of the evolution, the noise can safely be discarded. Thus stochasticity 
only manifests itself in the initial stages of the evolution, but in a crucial way. 
Indeed, as one can see in Fig. 15, after averaging, (nt) differs significantly 
from the mean-field result, and this difference stems from rare realizations in 
which the particle number stays low for a long time. Therefore, in individual 
realizations, stochasticity should accurately be taken into account as long as 
n t < n. Fortunately, when the number of particles in the system is small 
compared to the parameter N that fixes the typical maximum number of 
particles in a realization, the stochastic evolution is essentially governed by a 
linear equation. 
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Thanks to this discussion, we may assume that the evolution is linear as long 
as there are less than n particles in the system and deterministic when n t > n. 
It is then enough to compute Pn(i) for an evolution without recombinations, 
and (n t \n, t) for an evolution without noise. The second quantity is most eas- 
ily computed by replacing n t in Eq. (32) by the average quantity n^ F (or 
equivalently by discarding the noise term in Eq. (64)) and neglecting the term 
nf [F /A r (which is small compared to the term n^ F ). One gets a closed equation 



for n^ F in the form 



j MF (n MF ) 
u " / t MF \' l t ) 



(81) 



dt 1 N 
which is solved by 

= TT £^ = (n t \n,t) (82) 



n 



where the initial condition at time t = t has been chosen in such a way that 
nt = n. 

As for the distribution pnif) for the waiting times t to observe n particles in the 
system, its derivation is a bit more subtle. Since the evolution is taken linear 
until there are n particles in the system, the number of particles increases with 
time in any given realization. Then the following relation is true 



.-. d 
P ^ = dt 



E P M) (83) 



t=t n=n 



where P(n, t) solves the master equation (31) in which terms of order 1/N are 
discarded. This relation only holds because the probability that n be larger 
than n reads 

oo 

Prob(n >n,t) = J2 p ( n , (84) 

n=n 

thanks to the fact that n never decreases in realizations when nonlinear effects 
are neglected. 

We could solve the simplified equation for P(n, t), but for the sake of presenting 
a method that may be more general, we shall follow a slightly different route 
and establish first an equation that gives pn(t) more directly. 

Let us introduce Q(n,t) the probability that the number of particles remain 
strictly less than n for any time in [0,t], starting with a system of n particles 
at time 0. Then we obviously have 

jf di Pn (t) = Q(l,t). (85) 

which by simple derivation of Q(l, t) with respect to t gives the relevant distri- 
bution. We now establish an evolution equation for Q. Recall that the evolution 
equation for P was obtained by considering the variation in the number of par- 
ticles in the system between times t and t + dt. Here we consider the beginning 
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of the time evolution, between times and dt. The probability that the system 
does not exceed n particles up to time t + dt starting with n particles at time 
t = 0, Q(n,t + dt), is the probability ndt that the system gains a particle 
between times and dt multiplied by Q(n + l,t), plus a unitarity-preserving 
term. In this way, after having taken the limit dt — > 0, we get 

*?M = n(Q(n + M)-Q(M))- (86) 

This equation is valid when we neglect recombination processes, which is the 
relevant approximation here. In order to find a solution, we introduce the 
generating function for the moments of n: 

oo 

G(u,t) = ]>> n Q(n,t). (87) 

n=0 

The evolution of Q implies 

£ = (i-«)£-ia (88) 

at ou u 

This equation may be solved by the method of characteristics well-known 
for example in fluid mechanics, but also in QCD where it is commonly used 
to solve the renormalization group equation. We provide all details of the 
derivation of the solution in our simple case since it is not used so often in the 
particular field of high-energy QCD. 

The method consists in promoting the independent variable u to a function 
of time: u — > u(t). One then writes the total time derivative of G as 

dG(u(t),t) _ dG{u{t),t) dujt) dG(u(t),t) 

dt ~ dt + ~dt Ihi ' ^ ' 

Identifying the right handside of this equation to Eq. (88), one gets 

dGWM = _I G , (90 ) 
at u 

provided that u(t) solves 

d ^ = u-l. (91) 
dt v 1 

This equation is easily integrated: 

u(t) = 1 + (uo - l)e*, (92) 

where the initial condition u = u(0) is taken at zero time. The backward 
solution is also needed: 

u = 1 + (u(t) - l)e~*. (93) 
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Next, one integrates the ordinary differential equation (90) 

G(u(t),t) = G(u ,0)exp (^-j^dt'-^j (94) 

Replacing u(t') by its value given by Eq. (92) under the integration sign, then 
replacing u by its expression as a function of u and of t (Eq. (93)) one gets 

G(u,t) = G(l + (« - l)e-',0) -—. (95) 

U — 1 + 6 

Finally, the initial condition for G stems from the fact that Q(n,0) = for 
n > n and Q(n, 0) = 1 for n < n. Therefore, 

n-l i _ n 

G(u, 0) = 5> B = . . (96) 

n=0 1 - U 

Inserting this result into Eq. (95), we get 

G(M) » IziizOziOO!. 

1 — u 1 — e *(1 — -u) 
(5(1, t) is easily obtained from G, by a simple integration: 

where the integration runs over an appropriate contour in the complex w plane. 
We get from the Cauchy theorem 

Q(M) = 1 + (l - e-'X' 1 ' (99) 

and 

Pn(t) = = (#i _ i) e -f(i _ e -*f - 2 . (100) 

In the limits n ^> 1 and f ^> 1 which are relevant here, the distribution 
simplifies to 

p n (t) ~ ne _f - fie "'. (101) 
This is a Gumbel distribution. 

Plugging Eqs. (101) and (82) into Eq. (80), we get for the average number of 
particles after t units of time of evolution: 

roc _ frQ~t~™ e * 

(nt) = Nj Q dt i + ^_ {t _ iy (102) 



n 



Because the Gumbel distribution is strongly damped for t < 0, the lower 
integration boundary may safely be extended to — oo. Indeed, it is easy to see 
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that a conservative upper bound for the contribution of the domain ] — oo, 0] 
to the integral is e _n , which is very small in the limit n 1. Finally, we 
perform the change of variable b = ne~ l jr to arrive at the form 



It can be checked that it is exactly the form found through the diagrammatic 
approach to Pomeron field theory (compare Eq. (103) to Eq. (63)). 

The factorization in Eq. (80) and the convenient approximations that it sub- 
sequently allows are actually very important. Indeed, we realized that we may 
write the average number of particles at time t, whose expression would a pri- 
ori be given by the solution of a nonlinear stochastic differential equation, by 
solving two much simpler problems. The key observation was the following. 
When the number of particles in the system is low compared to the maxi- 
mum average number of particles N allowed by the reaction process, then the 
nonlinearity is not important, but the noise term is instead crucial. On the 
other hand, when the number of particles is large compared to 1, then the 
noise may be discarded, but the nonlinearity of the evolution equation, which 
corresponds to recombinations of particles, must be treated accurately. From 
this method, one gets an expression for (n t ) up to relative corrections of order 



When we address the problem of reaction-diffusion with one spatial dimension, 
we will rely on the very same observation. It is essentially the latter which will 
enable us to find analytical results also in that case. 

3.4 Relation to high energy scattering and the parton model approach 

So far, we have focussed on the factorial moments of the number n of particles 
in the system. We have seen how they may be computed from "Pomeron" 
diagrams, which are quite similar to the diagrams that appear in effective 
formulations of high-energy QCD. However, the relation to scattering ampli- 
tudes, which are the observables in QCD, may not be clear to the reader at 
this stage. In particular, we do not understand yet what would correspond to 
boost invariance of the QCD amplitudes. The aim of this section is to clarify 
these points. 

Let us consider a realization of the system of particles, evolved up to time t (at 
which it contains n t particles), that we may call the projectile. A convenient 
formalism to compute the weights of Fock states was presented in Sec. 3.2.1. 
We imagine that at time t, it scatters off a target consisting of a single particle, 
and can have at most one exchange with the target, which "costs" a factor 




(103) 



1/N. 
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1/N. All the particles in the system have an equal probability to scatter. Hence 
the probability that the system scatters reads T = n t /N. The average of T 
over events is the average particle number normalized to N. 

This way of viewing the evolution of the system makes it obviously very similar 
to the QCD dipole model introduced in Sec. 2, provided one identifies the 
number of particles to the number of dipoles and the time to the rapidity 
variable. The average of T over realizations is the elastic scattering amplitude. 

From this analogy, there is a property similar to boost invariance that should 
hold. Instead of putting all the evolution in the projectile, we may share it 
between the projectile and the target. Let us call n t i the number of particles in 
the projectile at the time of the interaction, and m t _ t i the number of particles 
in the target. The total evolution time is the same as before. To establish 
the expression for T in this frame, it is easier to work with the probability 
S — 1 — T that there is no interaction. If any number of interactions were 
allowed between each pair of particles from the projectile and the target, 
then one would simply write S = exp(— n t im t _ t i / N). But since the number of 
interactions should be limited to one per particle, one has to decrease n and 
m for each new power of 1/N, i.e. for each additional rescattering: 

1 11 

S = l- —nm + — — [ n (n ~ l)][m(m - 1)] 

- ^JpW n - " 2)]M™ " " 2)] • • • (104) 

where the time dependences are understood in order to help the reading. This 
is like a "normal ordering" of the expression to which we would arrive by 
assuming any number of exchanges. Note that S is not necessarily positive in 
a given event, and hence one looses the probabilistic interpretation once one 
has performed the normal ordering. 

Taking the average over realizations, one gets 

w = £((^),(<s^L*i^"- (105) 

If t' = t — t', the first two factors in each term of the series are of course 
identical after averaging. The sum runs over the number of actual exchanges 
between the probe and the target. A realization of the evolution, which would 
correspond to an event in QCD, is represented in Fig. 16. Note that the figure 
is very similar to Fig. 7a, except that particle mergings are allowed, while they 
have not been properly formulated in QCD yet. 

Now this expression should be independent of t'. It is not difficult to check that 
this is indeed true by taking the derivative of (S) with respect to t'. Expressing 
the averages of the factorial moments of the number of particles with the help 
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Fig. 16. Representation of the scattering of two systems of particles. The systems 
evolve in time from the left to the right. The horizontal lines represent the particles, 
and the vertical dashed lines the interactions between the systems. Each of the 
elementary scatterings comes with a power of 1/N. Note the strong similarity with 
the QCD diagram in Fig. 7a, except that in the present case, recombinations are 
included in the evolution of each of the systems. 

of the probability distributions P(n, t') and P(m,t — t') respectively, each term 
of the sum over k and m, n reads 



d(S) 



df 



= (PP -PP ) — m! W (106) 

, f. , [ n m n m) {n-k)\{m-k)\k\N k - 1 ' 

n,m,k fixed v ' y ' 



The time dependence is understood, and we introduced the notation P„ = 
P(n, •) and P n = dtP(n, •) to get a more compact expression. The time vari- 
abe that should be used for each factor is unambiguous since it is in one-to-one 
correspondence with the particle number index. We may use the master equa- 
tion (31) to express the time derivatives: 



p p —p p 



P m -\n <-> ml. 



(107) 

Recalling that there are sums over m, n and k which go from to oo, one 
may shift first the indices m and n in order to factorize P n P m in each term. 
The factors 1/N may then be absorbed by shifting k for the relevant terms. 
Then cancellations occur between the terms of both squared brackets in such 
a way that once the summations over n, m and k have been performed, the 
global result is 0. This proves the independence of (S) upon t', that is, "boost 
invariance" in a relativistic quantum field theory language. Of course, boost 
invariance is a consequence of some subtle interplay between the form of the 
interaction (104) and the form of the evolution encoded in the master equa- 
tion (31). Had we not normal ordered the expression for S in Eq. (104), boost 
invariance would not have hold as we shall check shortly. 
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We have seen that we may formulate scattering amplitudes in the zero-dimensional 
toy model, exactly in the same way as in QCD. We have seen in particular 
how crucial it is to include particle mergings consistently with the form of the 
interaction between the states of the projectile and of the target at the time 
of the interaction, in order to get a boost-invariant amplitude. 

3.5 Alternative models in dimensions 

For the sake of completeness, we shall now construct some variants of the 
zero-dimensional model introduced above, since the latter were also discussed 
in the literature. We review two of the most popular models. 

3.5.1 Allowing for multiple scatterings between pairs of particles 

Instead of assuming that there is at most one single exchange between each 
pairs of partons, one may allow for any number of exchanges. Then the defi- 
nition of S is modified as follows: 

(S) = (e~ s F)= Yl P(n,t')P(m,t-t')e-^. (108) 

n,m>l 

One sees immediately that if the probabilities P satisfy the master equa- 
tion (31), then this expression cannot be boost-invariant (i.e. independent of 
U). Indeed, if Eq. (31) holds, then 

P(n, t -> oo) = 5 n)N and P(n, t = 0) = 5 n>l . (109) 

It follows that in the frame in which the projectile is at rest, 

(S) t >= ^oo = e- 1 (110) 

while in the center-of-mass frame (if the projectile and the target share an 
equal fraction of the evolution), 

(S)*=l++oo = e- N (HI) 

which is very different. Actually, in this model, the average number of particles 
cannot saturate to a fixed value N. It would not be compatible with boost 
invariance. 

In order to preserve boost-invariance, one has to modify the master equation. 
We may write most generally 

Pn = T,(°LkPn-k-a£P n ). (112) 
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The coefficients a* are the transition rates from a (n — fc)-particle state to a 
n-particle state. We determine the a\ from the boost-invariance requirement. 
Actually, only one coefficient a^ =1 is needed in the case of this model. 

Using the same method as the one employed for checking the boost invariance 
in the previous model, we write 

f = Efe-^)(^). (H3) 

and express P n , P m with the help of the master equation (112). Requiring that 
the sum over n and m vanishes leads to the rates 

*J = N (l - e~ n ' N ) , (114) 

where the overall constant is determined from the rate in the unsaturated 
version of the model, which should hold for values of n <C N. This model was 
first proposed by Mueller and Salam [54]. 

We see that the saturation mechanism is quite different than in the previous 
model. Indeed, the average number of particles in the system keeps growing, 
but at a rate that slows down and depends on the number of particles in 
the system itself. Unitarity of the scattering probability T is ensured first by 
multiple scatterings rather than by the saturation of the number of particles 
to a constant number N (up to statistical fluctuations). 

This model was studied in detail in Ref. [95]. The conclusions drawn in there 
is that the saturation mechanism implied by the above model is likely to be 
quite close to the one at work in QCD. We could get analytical results for this 
model using one of the methods presented above. In particular, the statistical 
method outlined in Sec. 3.3 would apply and lead in a straightforward way to 
the expression for (n), up to corrections of relative order 1/N. 



3.5.2 Reggeon field theory 

Starting from the field theory formulation in Sec. 3.2.2, we may discard the 
4-Pomeron vertex (term (tf) 2 b 2 /N in Eq. (51)). The new Hamiltonian then 
reads 

H RFT = _ 6 f 6 _ ( 6 t)2 ft + J_ 6 t^ (n5) 

The stochastic formulation reads 

d z 

-^ = z t -^ + V^~ t vt + dt (116) 

(Compare to Eq. (73).) This is the zero-dimensional version of the stochas- 
tic equation defining the so-called Reggeon field theory, which was intensely 
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studied in the 70's as a pre-QCD model for hadronic interactions. 



This model has peculiar properties if one insists on interpreting it as a particle 
model. Indeed, the Hamiltonian (115) corresponds to a generating function for 
the factorial moments of the number n of particles in the system at a given 
time t that satisfies the partial differential equation 

dZ(z,t) _ ^ dZ(z,t) _ z &Z{z,t) 

dt ~ z ^ + z > dz N dz 2 ^ lU > 

and the corresponding master equation, obeyed by the probability P(n, t) to 
find n particles in the system at time t, writes 

^ = -nP(n, t) + (n- l)P(n - 1, t) 

+ l(n + l)(n + 2)P(n + 2, t) - ^n{n + l)P(n + l,t). (118) 

One can read off this equation the rates for particle creation/disappearance. 
One has a 1 — > 2 splitting, with rate dt; a 2 — » annihilation with rate dt/N; 
and a 2 — > 1 recombination with rate —dt/N. This is a negative number, and 
of course, it is unacceptable for a physical probability not to take its values 
between and 1. But we should not reject a priori negative probabilities as a 
formal calculation tool, as long as the physical probabilities are well-defined. 
However, a Monte-Carlo code based on these negative rates turns out to be 
extremely unstable, and thus of no practical use. 

Note that the statistical approach teaches us that in the iV 3> 1 limit, the 
moments of the number of particles in the system should not be very different 
than for the model with 3 and 4-Pomeron vertices, since it is essentially the 
form of the fluctuations in the dilute regime that determine the moments at 
all times. 

A detailed study of the special properties of this model as well as a comparison 
with reaction-diffusion-like models may be found in Ref. 



4 Review of general results on stochastic traveling-wave equations 



In Sec. 2, we have shown the relevance of the stochastic FKPP equation for 
high-energy QCD. The latter represents (classical) particle models that un- 
dergo a branching-diffusion process in one dimension, supplemented by a sat- 
uration mechanism. Sec. 3 was dedicated to a detailed study, from different 
points of view, of simplified models obtained from the former ones by switch- 
ing off diffusion. We now go back to the study of one-dimensional models. We 
proceed by steps: First, we shall address the deterministic FKPP equation 
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(which is equivalent to the BK equation in QCD) (Sec. 4.1). Second, we shall 
introduce fluctuations to get solutions for equations in the universality class 
of the sFKPP equation (Sec. 4.2 and 4.3). 

4-1 Deterministic case: the FKPP equation 

We address the simplest reaction-diffusion equation, namely the FKPP equa- 
tion 

d t u = d 2 x u + u - u 2 . (119) 

This equation was found to describe scattering in QCD under some assump- 
tions, see Sec. 2. 

It is a mathematical theorem [97] that this equation admits traveling waves 
as solutions, that is to say, soliton-like solutions such that 

u(t,x) =u(x-vt) (120) 

where v is the velocity of the wave, u is a front that smoothly connects 1 
(for x — > — oo) to (for x +oo). The velocities of the traveling waves and 
their shapes for large x are also known mathematically. Starting with some 
given initial condition which itself is not necessarily a traveling wave such as 
Eq. (120), the solution converges at large times to a stationary wave front. 
The front velocity during this phase may also be predicted asymptotically. 
We informally review these results in this section. 

4-1.1 General analysis and wave velocity 

The FKPP equation (119) encodes a diffusion in space (through the term d 2 u 
in the right handside), a growth (term u), and a saturation of this growth 
(term — u 2 ). It admits two fixed points: the constant functions u(t,x) = 
and u(t,x) = 1. A linear stability analysis shows that is unstable, while 1 
is stable. Indeed, thanks to the growth term u in the right handside, a small 
perturbation u(t,x) = £ C 1 grows exponentially with time. On the other 
hand, a perturbation near 1 of the form u(t, x) = 1 — e goes back to the fixed 
point 1 through evolution. Hence the FKPP equation describes the transition 
from an unstable to a stable state. Therefore, we expect that the linear part 
of the equation drives the motion of the traveling wave, since the role of the 
nonlinear term is just to stabilize the fixed point u—1. 

We shall cast the linear part of the equation in a more general form: 

d t u(t,x) =u(-d x )u(t,x), (121) 



49 



where u(—d x ) is a branching diffusion kernel. It may be an integral or differ- 
ential operator. An appropriate kernel is, in practice, an operator such that 
the "phase velocity" v (7) —oj{ r y)/ r ) (see below) has a minimum in its domain 
of analyticity. The FKPP equation corresponds to the choice lo(— d x ) = d% + l. 

Let us follow the wave front in the vicinity of a specific value of u. To this 
aim, we define a new coordinate xwf such that 

x = %f + vt. (122) 

The solution of the linearized equation (121) writes most generally 

f ^7 

u(t, x) = J c — : uo(t) ex P (-7(^wf + vt) + w(i)t) , (123) 

where a; (7) is the Mellin transform of the linear kernel cj(— d x ) (and thus 
7 corresponds to —d x ), and defines the dispersion relation of the linearized 
equation. 1x0(7) i s the Mellin transform of the initial condition u(t = 0,x). 
Let us assume that the initial condition is a function smoothly connecting 
1 at x — —00 to at x = +00, with asymptotic decay of the form u(t = 
0,x) ~ e~ loX . Then 1x0(7) nas singularities on the real negative axis, and on 
the positive axis starting from 7 = 70 and extending towards +00. Let us 
take a concrete example: If u(0,x < 0) = 1 and u(0, x > 0) = e^ lox , then 
1x0(7) = 1/7 + l/(7o — 7)- The integration contour C should go parallel to the 
imaginary axis in the complex 7-plane and cross the interval [0, 70]. 

Each partial wave of wave number 7 has a phase velocity 

vM = ^, (124) 

7 

whose expression is found by imposing that the exponential factor in the 
integrand of Eq. (123) be time-independent for v = 1^(7). 

We are interested in the large-time behavior of u(t,x). The integrand in 
Eq. (123) admits a saddle point at a value 7 C of the integration variable such 
that 

w'(7c) = v, (125) 
that is to say, when v coincides with the group velocity of the wave packet. 
But the large-time solution is not necessarily given by the saddle point: This 
depends on the initial condition 1x0(7). 111 order to understand this point, let 
us work out in detail the simple example of initial condition quoted above. 
The integral has two contributions for large t: 

u (j- £j — e -7o(a^VF+ft)+w(7o)* _|_ Ke -lc(x W F+vt)+uj("/ c )t (126) 

up to a relative 0(1) factor k. The time invariance of u(t,x) in the frame of 
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the wave may only be achieved by tuning v to one of the following two values: 



(ii) v c 



^(7o) 

7o 
w(7c) 



w'(7c) 



(127) 



In the second case, i> coincides with the minimum of the phase velocity oj{ r y)/ r ) 
and in particular, v c < Vq. The relevant value of v depends on the shape of 
the initial condition: 

• If 70 < 7c, i.e. the decay of the initial condition is less steep than the decay 
of the wave from the saddle-point, then one has to pick the first choice (i) for 
the velocity. Indeed, this is the only one for which the first term in Eq. (126) 
is time-independent, and the second term vanishes at large time. Due to the 
fact that v c < vq, choice (ii) would make the first term in Eq. (126) blow 
up exponentially, u ~ e y °( v °~ Vc ^ t . 

• If instead 70 > 7 C , then it is the second choice (ii) that has to be made. The 
saddle point dominates, and the wave velocity at large time is independent 
of the initial condition. 

Fig. 17 summarizes these two cases. 

The limiting case 70 = 7 C requires a special treatment. Since it is not relevant 
for the physics of QCD traveling waves (only the case 70 > 7 C is actually 
relevant), we refer the interested reader to the review paper of Ref. [72] for a 
complete treatment also of that case. 

There exists a rigorous mathematical proof of these solutions in the case of the 
straight FKPP equation [97] . These results are largely confirmed in numerical 
simulations for various other branching diffusion kernels, including the ones of 
interest for QCD (see e.g. [98,74], and Ref. [99,33,100] for earlier simulations 
of the BK equation). 

Actually, in QCD as well as in many problems in statistical physics, the initial 
condition is localized or has a finite support, and hence, its large- a; decay is 
always very steep. Thus for the physical processes of interest in this review, 
the asymptotic front velocity, that we will denote by for reasons that will 
become clear later, reads 



where the last equality defines 7 C . Note that in the context of particle physics, 
this result was already known from the work of Gribov, Levin, Ryskin [13], 
and was rederived later in the framework of the BK equation [34,101,35]. 



V OQ =v c = 



^(7c) 
7c 



(128) 
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Fig. 17. Front velocity as a function of its asymptotic decay rate 70 (dashed curve). 
It has a minimum at 7 = j c . The full line represents the actual velocity that 
would be selected starting with an initial condition decaying as e~ lox for large x. If 
70 = 7_ < 7 C (initial condition less steep than 7 C ), then the asymptotic velocity is 
the phase velocity of a front which has the same asymptotics as the initial condition. 
For any 70 = 7+ > 7 C , the velocity of the front is the minimum of the phase velocity 

So far, we have discussed the asymptotic velocity of the solutions to the FKPP 
equation as a function of the initial condition. When the initial condition is 
steep enough, then the asymptotic front velocity takes a fixed value which is 
the minimum of uo( r ))/ r y. In the opposite case, the shape of the initial condition 
is retained (see Fig. 18). We wish to know more detailed properties of the wave 
front, such as its shape and the way its velocity approaches the asymptotic 
velocity. There are several methods to arrive at this result. At the level of 
principle, they all rely on a matching between a solution near the fixed point 
u—1, and a solution of the linearized equation which holds in the tail m C 1. 



4-1.2 Diffusion equation with a boundary 

We now come back to the original FKPP equation (119). We have seen that 
the nonlinearity — u 2 has the effect of taming the growth induced by the linear 
term u, when u gets close to 1. But nonlinear partial differential equations are 
very difficult to address mathematically. It may be much simpler to address 
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Fig. 18. Sketch of the shape of the front according to the large- x behavior of the 
initial condition u(t = 0, x) ~ e~ li,x . Top: 70 < 7 C . The asymptotic shape of the ini- 
tial condition is conserved. The relaxation of the front is fast. Bottom: 70 > j c - The 
asymptotic shape of the front is e~^ cX , and the velocity for t = 00 is v c = uj(^ c )/^ c . 
The asymptotic shape is reached over a distance \/t ahead of the front, and the 
velocity at finite time is less than the asymptotic velocity by 277- 

the linear equation 

d t u = d 2 x u + u (129) 
supplemented with an absorptive (moving with time) boundary condition that 
ensures that u(t,x) has a maximum value of 1 at any time. We need to work 
out the solution of Eq. (129) with this kind of boundary condition. Here, 
we reformulate the approach proposed in the QCD context by Mueller and 
Triantafyllopoulos [35] (see also Ref. [102] for an account of the next-to-leading 
order BFKL kernel). 

A solution to Eq. (129) with initial condition u(t — 0, x) — 5(x — x ) is given, 
for positive times, by 

^^vs^t'-^)' (130) 

This solution holds if the boundary condition is at spatial infinity. The solution 
of the pure diffusion equation, without the growth term, is of course nothing 
but u(t,x)e~ t . We shall denote it by -u PD (t,a:). 

If instead of the boundary condition at infinity there is an absorptive barrier 
at say x — X, i.e. if u(t,x = X) = for any t, then a solution may be 
found through a linear combination of the latter solution with different initial 
conditions, in such a way as the sum vanishes at x — X. This is known as the 
method of images. It is based on the observation that any linear combination 
of Eq. (130) also solves Eq. (129). From the solution with initial condition 
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Fig. 19. Shape of the solution of the branching diffusion equation (129) with a 
moving cutoff, whose position is adjusted in such a way that the maximum of 
u(t, x) be 1 at all times. The solution is represented at two different times t\ and 
t2, showing the soliton-like behavior of the solution. 

S(x — x ), we subtract the solution of the same equation but with initial 
condition 5(x — (2X — xq)) , in such a way that the solution vanishes for x = X, 
at any time. We get 

e* / (x-xpf (x-2X+xq) 2 

ux(t, x) = ■ e « — e 4* 

We do not expect the solution to this problem to represent accurately the 
solution to the full FKPP equation near the boundary x ~ X. So the region 
of interest will be ahead of the boundary by a few units, while the starting 
point Xq of the evolution is at some finite distance of the boundary: 

x-X>l and x - X ~ 1. (132) 

One may then expand the two Gaussian terms: 

x -X x-X f {x-Xf \ 

Uxi ,x) = exp y — • (133) 

The solution to the simple diffusion equation without the growth term, namely 
d t u = d^u, is the one that we will actually use in the following. It would again 
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be the latter solution scaled by e *, namely 

x — X x — X 



u™(t,x) 



An~ * 3 / 2 



cxp 



[x-Xf 
At 



(134) 



where the superscript PD stands for "pure diffusion" . Note that in this equa- 
tion, X does not depend on time. We cannot implement in a straightforward 
way a time-dependent absorptive boundary. We will get to such a solution by 
successive iterations: The main trick is to go to a frame in which the solution 
of the branching diffusion with a boundary is stationary for large times. 

Let us start from the solution uq in Eq. (130). The lines x of constant uo(t, x) = 
C (without a boundary) are obviously given by 



x = x + It — - \ogt — log(CV47r) + terms vanishing for t — > oo. (135) 

(We have selected the rightmost solution x > x ). Let us change frame by 
writing x = x\ + xq + 2t. Then in this new variable, u(t, x) in Eq. (130) reads 



u(t, x) — e 



-X! 



e 4t 



/ 47rt 



= e - Xl u PD (t,a;i), 



(136) 



where we have factored out the solution of the pure diffusion equation, but this 
time, in the moving frame defined by the coordinate x±. We may implement 
an absorptive boundary condition, fixed in this new frame, by replacing -u PD 
by w PD in Eq. (134). Note that X is fixed with respect to x 2 , but in the 
original frame defined by coordinate x, it is of order 2t. The solution has lines 
of constant u which solve 



Xi = --log* 



^ + log (xi - X) - log —t— . 

x - X 



At 



(137) 



The two last terms are subdominant because according to Eq. (135), X\ — X 
log*, and because xq — X is a constant. We further define a new frame 



3 , , CVAn 
x 2 = x 1 + - log* - log 

2 xq — X 



X 



Xq 



2*+ -log* 



l ° g x^X (138) 



Going back to the expression for u (t,x) (see Eq. (130)), we substitute x by 
its expression as a function of x 2 and get 



u 



(t,x) =t 3/2 e~ X2 



CVAn 
x -X 



e « 



'Ant 



x exp | — — h x 2 



At 



9 log* + log^ , 
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We replace the expression inside the squared brackets, which is nothing but 
■u PD , by the solution with a boundary m™. We check that the value of x 2 for 
which u(t,x 2 ) is constant is now a mere constant for large t. Going back to 
the original frame, we get 

u(t, x) = Ce- X (x - X(t))e-( x - x ^ exp (^- {x - f^ ^j , (140) 

where 

X(t) = 2t-^logt + 0(l) (141) 

is the position of the absorptive boundary for large times, and thus, the po- 
sition of the front. The constant X is the position of the front in the moving 
frame. Setting X = — 1 and C — 1, the maximum of u is reached at x — X{t), 
and is indeed equal to 1. 

For large t or in the region x — X(t) < \fi which expands with time, the 
Gaussian factor goes to 1, and we see that u(t,x) only depends on one single 
variable x — X(t). This was expected: It is precisely the defining property of 
traveling waves. But in addition to these asymptotic solutions, we get from this 
calculation the first finite-t correction to the front shape and front velocity. 

Actually, the speed of the front is intimately related to its shape. At time t, 
it has reached its asymptotic shape over the distance y/i from the saturation 
point. This remark will be important in the following. 

We have derived the solution of a problem that was not exactly the initial one, 
however, we believe that the shape of the front in its forward part (u< 1) 
as well as its velocity are quite universal. Indeed, physically, these properties 
are completely derived from the linear part of the equation. For this reason, 
the front is said to be "pulled" by its tail. The nonlinearity only tames the 
growth of u near u ~ 1, and so its precise form should not influence the 
front position itself, at least at large enough times. Thus we expect these 
solutions to have a broad validity, only depending on the diffusion kernel, 
and so, may be obtainable from our calculation up to the replacement of 
the relevant parameters. For the more general branching diffusion kernel in 
Eq. (121), the velocity of the front would read 



dX 


_ w(7c) 


3 


~dt 


7c 


2 7c t + 



(142) 



where 7 C solves o;(7 c ) = Jc^'ilc), as was explained in Sec. 4.1.1. The front 
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shape in its forward part x — X(t) ^> 1 reads 



u(t, x) = (x- X(t))e'^ [x - X(t)) exp 



2u»{ lc )t 



(143) 



up to an overall constant. Fig. 19 represents a sketch of the solution at two 
different times. Note that the asymptotic shape is an exponential decay, 

u(t,x) ~ e -^ {x - x{t)) . (144) 

From Eq. (143), this shape extends over a range 



L — x — X(t) ~ J2u"(-f c )t. (145) 



In other words, the time needed for the front to reach its asymptotic shape 
over a range L reads 

L 2 

t — . (146) 



Through our simple calculation, we got the lowest order in an expansion of the 
front shape and position at large times. The next corrections to X(t) would 
be of order 1 (this constant depends on the way we define the position of the 
front), followed by an algebraic series in t whose terms all vanish at large t. The 
first next-to-leading term in the series has been computed (see Ref. [103]): It 
turns out to be of order l/y/i. We will not reproduce the calculations that lead 
to it because they are rather technical and there is already a comprehensive 
review paper available on the topic [72]. But let us write the result for the 
position and the shape of the front at that level of accuracy, for the more 
general branching diffusion kernel given by Eq. (121). To that accuracy, the 
front position reads [103] 



For the simple FKPP case, we recall that u(— d x ) = d 2 + 1, then 7 C = 1 and 
^(Tc) — 2. The first two terms in the last equations match the ones found in 
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Eq. (142). The shape of the front in its forward part has the following form: 

i(t,x) = C ie -^ X ~ X W cxp (-z 2 ) 



X 



7c (x - X{t))\ + C 2 + 3 - 2C 2 + 



w"(7c) , 



'2 7c" (3) (7c) 
.3 w"( 7c ) ' 3 

+ 6^ (l - 1F1 



I? lj §5 3; ^ 



1 3. r ^ 

2)2' 



where 



z = 



x - X{t) 



)z + 0(l/Vt)\, (148) 



(149) 



and 2F2, 1F1 are generalized hypergeometric functions. The terms in the first 
line match with the result of our calculation (Eq. (140)) for the relevant value 
of 7 C . These expressions should apply also to QCD, up to the relevant replace- 
ments given in Tab. 1. 

So far, we have considered equations of the type of Eq. (119) as saturation 
equations, in the sense that they describe the diffusive growth of a continuous 
function u until it is tamed for u ~ 1. We will see below that these equations 
may actually be given a different physical interpretation. 



4-1.3 Discrete branching diffusion 

We have investigated the solutions of the FKPP equation in a mathematical 
way, without discussing the physics that may lead to such an equation. The 
absorptive boundary that we have put replaces the nonlinear term in the 
FKPP equation, whose role is to make sure that u never exceeds the limit 
u — 1. Hence we have thought of this boundary as a way to enforce the 
saturation of some density of particles. Actually, the FKPP equation (119) may 
stem from a branching diffusion process in which the number of particles is 
unlimited, and thus, for which there is no saturation at all. As a matter of fact, 
this is the way how the BK equation is built in QCD: An exponentially growing 
number of dipoles, stemming from the rapidity evolution of a hadronic probe, 
scatters off some target. The overall interaction probability is unitary because 
multiple scatterings are allowed (the interaction probability of n dipoles is 
actually of the form e - " 3 ™), but not because there is a saturation of the number 
of dipoles in the wavefunction of the probe. We refer the reader to Fig. 5 for 
a picture of the process. 

To illustrate how the FKPP equation arises in such a simple model of branch- 
ing diffusion, let us consider a set of particles on a line, each of them being 
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Fig. 20. Example of branching diffusion process on a line (see the text for a mathe- 
matical description of the evolution rules). If the number of individuals is limited by 
a selection process which, at each new branching, eliminates the individual sitting 
at the smallest x as soon as the total number of individuals reaches say N (N = 10 
in this figure), then only the branches drawn in thick line survive. 

indexed by a continuous variable x. (Such a model was considered for instance 
in Ref. [71]). We let the system evolve according to the following rules. During 
the time interval dt, each particle has a probability dt to split in 2 particles. 
Unless it splits, it moves of the small random amount 5x, which is a Gaussian 
random variable distributed like 



p(5x) 



1 



VAirdt 



exp 



(5x) 



Adt 



(150) 



Let us consider the number of particles n(t, x) contained in an interval of 
given size Ax centered around the coordinate x. At time t = 0, the system is 
supposed to consist in a single particle sitting at the origin x = 0. A sketch of 
a realization of this model is shown in Fig. 20. From the evolution rules, we 
easily get an equation for the average number of particles (n): 



(n{t + dt, x)) = dt 2(n) + (1 - dt) / d(5x)p(5x) (n{t, x - Sx)) 



:i51) 



which reads, after replacing p by Eq. (150) and after the limit dt — ► has 
been taken, 



9(n) d 2 (n) 
(n) + 



(152) 



dt dx 2 

All the dependence on the size Ax of the "bin" is contained in the initial 
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condition. It is clear that for large enough times, the solution to this equation 
is given by Eq. (130), and thus the lines of constant (n) are given by Eq. (135). 

Let us now define 

S(t,x) = e~ n ^ N (153) 

where N is some (large) constant. This definition is reminiscent of the S- 
function, related to the scattering amplitude, introduced in the discussion 
of the BK equation in Sec. 2. For large enough x, n(t,x) <C N and thus 
1 — S(t,x) ~ n(t,x)/N — > 0. For any x, the exponential makes sure that S 
ranges between and 1. Thus S (or 1 — S) has the shape of a traveling wave. 
Its position X(t) is the value of x for which n(t, x) is some given constant say 
of the order of N . Hence, up to fluctuations, it is given by Eq. (135). 

On the other hand however, the average of S over events, namely A = 1 — (S) 
obeys the FKPP equation. Indeed 

(S(t + dt,x)) = dt(S{t,x)) 2 + (1 -dt) f d(5x)p(5x)(S(t,x-5x)). (154) 

In the limit dt — > and rewriting the equation for A, we get 

OA _ cPA a _ a2 
dt dx 2 

Hence A is a traveling wave at large times, and its position X(t) is given by 
Eq. (141). It is obviously behind by a term logt with respect to the value of 
x for which the average number of particles has a given constant value. Fur- 
thermore, the probability distribution of the position of the rightmost particle 
(or of the k-th rightmost particle for any given k) may also be derived from 
the FKPP equation. It turns out that in any event, the average x for which 
n(t,x) has a given value, say no, moves with the FKPP velocity which can 
be read off from Eq. (141). This is much slower than the velocity with which 
X(t) defined in such a way that (n(t,X(t))) = n moves. 

All this may seem a bit paradoxical. But actually, it is just related to the fact 
that (e~ n / N ) cannot be approximated by e~( n ^ N . We may understand it in 
the following way. By taking the average of n, we have somewhat forgotten 
a fundamental property of n: its discreteness. Indeed, it only takes integer 
values, and in particular, the distribution of n in a realization has a finite 
support: At any time, there is a value of x to the right of which there are no 
particles at all. n obeys a stochastic equation. This is not the case for (n), 
which just obeys an ordinary branching diffusion equation. 

In order to recover the effect of the discreteness of n and compute the veloc- 
ity, we may again use the absorptive boundary trick. Let us solve the linear 
equation 

d t (n) = d 2 x (n) + (n) (156) 
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(155) 



<n> 




x 



Fig. 21. Solution of the branching diffusion equation (156) with a moving absorptive 
boundary that forces (n) to vanish at some well-chosen point. Two different times 
are represented. 

with an absorptive boundary. The absorptive boundary will be placed in such 
a way that at a distance of order one to its left (we will focus on the right- 
moving wave), (n) = 1 (see Fig. 21). There is no difference in principle with 
the boundary calculation that we have performed before, except that the ab- 
sorptive boundary is now placed to the right of the front (i.e. Xq < X in the 
notations used above). Thus we find without any further calculation that the 
realizations of n move, on the average, with the FKPP velocity (142). 

4-2 Combining saturation and discreteness 

We have seen that physically, the KPP equation (or the BK equation in QCD) 
may be interpreted either as an equation for the growth, diffusion and satura- 
tion of a continuous function, or as the evolution equation for the average of a 
bounded function of a discrete (thus stochastic) branching diffusion process. 
For each of these interpretations, we may find the main features of the solu- 
tions by imposing one absorptive boundary on the linear partial differential 
equation encoding branching diffusion. In one case, the boundary is a cutoff 
that prevents u to be larger than 1: It represents saturation, i.e. the explicit 
nonlinearity present in the FKPP equation. In the other case, the boundary 
forces the function n that represents the number of particles to vanish quickly 
when it becomes less than 1. Formally, it actually models the intrinsic discrete- 
ness of the number n of particles, and avoids to address a stochastic equation 



61 



fraction of particles to the right of x 



1/N 

x 

Fig. 22. Branching diffusion model of Sec. 4.1.3 with selection that limits the total 
number of particles to N. One sees that the fraction of particles to the right of x 
looks like a traveling wave front. 

directly. 

In physical cases such as reaction-diffusion processes for finite N, we define 
u(t,x) as the number of particles per site (or per bin) in x normalized to N. 
Hence it takes discrete values: 1/N, 2/N etc... While for large N discreteness 
is unlikely to play a role in the region u ~ 1, it is expected to be crucial 
when u ~ 1/N. It is thus natural to impose the two boundaries: one repre- 
senting saturation of the particle number, the other one discreteness of the 
same quantity. A model that these two cutoffs may represent is for example, 
the branching diffusion model in Sec. 4.1.3, but in which the total number of 
particles is limited to N by keeping only the N rightmost ones at each new 
branching. It is clear that the function U(t,x) defined to be the number of 
particles to the right of some position x normalized to the maximum number 
N is, for large enough times, a front connecting 1 (for x — > — oo) to (for 
x -> +oo) (see Fig. 22). 

React ion- diffusion problems (described by nonlinear stochastic partial differ- 
ential equations) were interpreted as branching diffusion problems taking place 
between two absorptive boundaries for the first time by Brunet and Derrida 
in Ref. [f 05] and, independently, by Mueller and Shoshi, in the case of QCD 
in Ref. [39]. Note however that the present interpretation of the cutoffs was 
only found in Ref. [40] in the context of the QCD parton model. Mueller and 
Shoshi considered both cutoffs for reasons related to the boost-invariance of 
the QCD amplitude. The duality of the two boundaries, that is to say of the 
dense and dilute regimes of the traveling wave, was studied more deeply in 
Refs. [106,107,108,109,110]. 

Before moving on to the technical derivation of the shape and position of the 
front in this case, let us figure out what we expect to find. 
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Starting from the initial condition, the front builds up and its velocity in- 
creases with t (see Eq. (142)) until it reaches its asymptotic shape, which is a 
decreasing exponential e~ Jc ^ x ~ x ^ that holds for all x ^> X(t). But if the front 
is made of discrete particles, then it has a finite support, and the exponential 
shape may not extend to infinity to the right, since u(t,x) has to be either 
larger than 1/N, or zero. It cannot take values that would be a fraction of 
1/N in realizations, and thus, we cannot accommodate the shape e~' 1c( - x ~ x ^ 
for arbitrarily large values of x, since it would mean authorizing arbitrarily 
small positive values of u(t,x). From Eq. (146) and from the shape of the 
asymptotic front (144), the exponential shape sets down to u = 1/N at time 

Beyond, the front cannot develop any longer, and thus, its shape and velocity 
remain fixed. t re iax is the time that is needed for the front to relax from any 
perturbation, which is why we have put the subscript "relax". 

From Eq. (142) evaluated at t — t rc i aX ) we get the new asymptotic velocity, 
that takes into account the effects of discreteness, in the form 

^ = ^M-3c 7 < (7c) . (158) 
dt 7c log N v ; 

The calculation of c requires a proper account of the exact shape of the front. 
We shall turn to this calculation now. 

As announced, we are now going to solve the linear branching diffusion equa- 
tion with two absorptive boundaries: one representing saturation, the other 
one discreteness. First, as in the one-boundary case, let us solve the simple 
diffusion equation d t u = d 2 x u between two boundaries, at X and Y respec- 
tively, that is to say, with the conditions u(t, X) = u(t,Y) = 0. The simplest 
method in this case is to take the ansatz 

u p x D Y (t,x) = f(t)g(x). (159) 

Then the diffusion equation reads 

/'(*) g"(x) 



fit) 9(x) 



= -A (160) 



A is necessarily a constant, being both a function of t only and of x only. The 
equations for / and for g are easily solved. All in all, we get for u 



u 



PD 

XY 



(t, x) = Ae' xt sin [V\(x - X )] , (161) 



where A and X Q are constants which we will shortly determine from the bound- 
ary condition. Note that only positive values of A are physical, since negative 
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u(t,x) 




X(t : ) X(t 2 ) 



X(tj)+L X(t 2 )+L 



Fig. 23. Solution to the branching diffusion equation with two boundaries. 

ones would lead to an exponential increase of the solutions. The boundary 
conditions at X and Y fix Xq to X and lead to a quantization of A: 



L 2 



A 



(162) 



where L = Y — X is the size of the wave front and k is an integer. The general 
solution is a sum of u^P Y over an possible values of k, with coefficients fixed 
by the initial condition. But at large time, thanks to the exponential decay of 
f(t) with t, only the mode k = 1 survives. The final solution thus reads 



= A exp I - 
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sm 



7l(x - X) 



(163) 



where the constant A is determined from the projection of the initial condition 
on the fundamental mode of the "cavity" [X, X + L] . 

We now need to determine the time dependence of X. It will follow from the 
search of the frame in which the front is stationary in time. 

The first step (determination of x±, see Eq. (136)) is the same as in the 
one-boundary case. Starting from Eq. (136), we substitute M PD (t,xi) with 
u^P x+L (t,xi) and look for the lines of constant u. This leads us to introduce 



7T 

X 2 = Xx + —t. 

Going back to the original variables, we get, for L large, 



(164) 



u(t, x) = Ae- X e~ {x - X{t)) sin £ {x - X(t)) 



(165) 
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where 



7T 



X(t) = 2t--t. 



(166) 



We are left with the determination of the size L of the front. Near the left 
boundary (at a distance of order 1), u should be of order 1, while close to the 
right boundary, it should approach 1/N. We write 



u(t,X(t) + l) = 1, u(t,X(t) + L-l) 



1 

N' 



(167) 



Then we see that L = \ogN and A = kL, where k = 0(1). The position of 
the rightmost boundary (Y) is the point to the right of which there are no 
particles in typical individual realizations. We will denote it by x tip (t). 

All in all, writing it for a more general diffusion equation d t u = cu(—d x )u, the 
final solution reads 



u(t,x) = Ke-^ x - x WLsm n{x ~X {t)) 



(168) 



(see Fig. 23) where the size of the front is 

logiV 



7c 



(169) 



and its velocity reads 



at 



7rV(7c) w(7c) tt 7c^"(7c) 



2 log 2 N 



(170) 



The subscript BD stands for "Brunet-Derrida" . For a; (7) = 7 2 + 1, 7 C = 1, 
a;(7c) = w"(7c) = 2 and we recover Eqs. (165), (166). 



^.5 Beyond the deterministic equations: Effect of the fluctuations 



So far, we have actually solved deterministic equations although we were ad- 
dressing a model with a discrete number of particles, that therefore has nec- 
essarily fluctuations. Our procedure gave the leading effects. We shall now 
incorporate more fluctuation effects, in a phenomenological way. (We shall 
essentially review Ref. [111]). 
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u(t,x) 




X 



L 




Fig. 24. Evolution of the front with a forward fluctuation. At time to, the primary 
front extends over a size L and is a solution of the branching diffusion equation 
with two appropriate boundaries. An extra particle is stochastically generated at a 
distance 5 with respect to the tip of the primary front. At a later time, the latter 
grows deterministically into a secondary front that is a bit slower, and that will add 
up to the primary one. The overall effect, after relaxation, is a shift to the right of 
the distance R(S) with respect to the position of the front if a fluctuation had not 
occured. 

4-3.1 Phenomenological model and analytical results 

The two-boundary procedure has led to the following result: The front prop- 
agates at a velocity Vbd in Eq. (170) lower than the velocity predicted by 
the mean-field equation (142), and its shape is the decreasing exponential 
e -ic(x-x(t)) down to the position 



(up to a global constant independent of N) at which it is sharply cut off by an 
absorptive boundary. This boundary was meant to make the front vanish over 
one unit in x, hence to implement discreteness on a deterministic equation. 

But since the evolution is not deterministic, it may happen that a few extra 
particles are sent stochastically ahead of the tip of the front (See Fig. 24). 
Their evolution would pull the front forward. To model this effect, we assume 
that the probability per unit time that there be a particle sent at a distance 
S ahead of the tip simply continues the asymptotic shape of the front, that is 
to say, the distribution of S is 



\ogN 



(171) 



p(5) = de"^, 



(172) 
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where C\ is a constant. Heuristic arguments to support this assumption were 
presented in Ref. [111]. Note that while the exponential shape is quite natural 
since it is the continuation of the deterministic solution (143) in the linear 
regime, the fact that C\ need to be strictly constant (and cannot be a slowly 
varying function of 5) is a priori more difficult to argue. 

Once a particle has been produced at position x t i P + 6, say at time to, it starts 
to multiply (see Fig. 24) and it eventually develops its own front (after a time 
trciax of the order of L 2 ), that will add up to the deterministic primary front 
made of the evolution of the bulk of the particles. 

Note that the philosophy of our phenomenological approach to the treatment 
of the fluctuations is identical to the spirit of the statistical approach in Sec. 3.3 
developped for the zero- dimensional model. Whenever the number of particles 
is larger than n (n — 1 here), we apply a deterministic nonlinear evolution. 
Fluctuations instead are produced with a probability which stems from a linear 
equation. 

Let us estimate the shift in the position of the front induced by these extra 
forward particles. Between the times to (of the order of 1) and t = t + t re i ax , 
the velocity of the secondary front is given by Eq. (142). Hence its position 
X^(t), after relaxation, will be given by 

X& (t) = x BD (t) + 6+ f dt' v t ,- t0 ~ X BD (t) + S-^- log L 2 (173) 

Jt 27c 

where Xbd(^) — Vbd^- Eq. (173) holds up to a constant independent of 5 and 
N. We have used Eq. (142) to express v t /- to . The observed front will eventually 
result in the sum of the primary and secondary fronts, after relaxation of the 
latter. Its position will be X BD (t) supplemented by a shift R(S) that may be 
computed by writing the resulting front shape in the large-x tail as the sum 
of the primary and secondary fronts: 

e -7c(x-X B D(t)-.R(5)) = e -7c(x-X BD (t)) _|_ e -lc(x-X{t)) 

where C2 is an undetermined constant. From Eq. (174) we get the shift 

1 / e 7c<5 \ 

R(6) = -log 1 + C 2 — . (175) 



The probability distribution (172) and the front shift (175) due to a forward 
fluctuation define an effective theory for the evolution of the position of the 
front X(t): 

X(t + dt) = l X{t) + VBDdt Pr ° ba - 1 ~ dt fo°° d5p(5) 

\X(t) + V BD dt + R(5) proba. p(S)dSdt. 



67 



From these rules, we may compute all cumulants of X(t), by writing the 
evolution of their generating function, deduced from the effective theory (176): 

JUog(e^) = A\/ BD + / dSp(S) (e XR ^ - l) (177) 

The left hand-side is a power series in A whose coefficients are the time deriva- 
tives of the cumulants of X(t). Identifying the powers of A in the left and right 
handside, we get 



V - Vbd = / d5p(5)R(5) 



dC 2 31ogL 

7c 7c^ 3 



[n-th cumulant] r C1C2 n\((n 



J d5p(5)[R(5)] n 



'178) 



t J 7 C 7™L 3 



We see that the statistics of the position of the front still depend on the 
product C1C2 of the undetermined constants C\ and Ci- We need a further 
assumption to fix its value. 

We go back to the expression for the correction to the mean-field front velocity, 
given in Eq. (170). From the expressions of R(S) (Eq. (175)) and of p(S) 
(Eq. (172)), we see that the integrand defining V — Vbd i n Eq- (178) is almost a 
constant function of S for 8 < 5 = 3 logL/7 c , and is decaying exponentially for 
5 > 5 . Furthermore, R(5 ) is of order 1, which means that when a fluctuation 
is sent out at a distance 8 ~ 5 ahead of the tip of the front, it evolves into a 
front that matches in position the deterministic primary front. We also notice 
that when a fluctuation has 5 < 5$, its evolution is completely linear until it 
is incorporated to the primary front, whereas fluctuations with 5 > <5 evolve 
nonlinearly but at the same time have a very suppressed probability. We are 
thus led to the natural conjecture that the average front velocity is given by 
Vbd in Eq. (170), with the replacement 

L^L cS = ho = h3 , (179) 

7c 7c 7c 

namely 

V = ^ ff . (180 ) 

7c ( 1okN + 31 °K 1 °K Jv r 

The large- N expansion of the new expression of the velocity yields a correction 
of the order of log log N/ log 3 N to the Brunet-Derrida result, more precisely 



T/ w(7 c ) tt 2 7 c w"(7c) . 2 «/ ^ 3 log log 

V = o h vr 7 c u; ( 7 C ) n . 

7c 2 log N ' KI J log N 



(181) 



Eqs. (178) and (181) match for the choice CiC 2 = n 2 uj"('y c ). From this de- 
termination of CiC 2 , we also get the full expression of the cumulants of the 
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position of the front: 



[n-th cumulant] 
t 



= vr 2 7 >"( 7c 



n!C(n) 
7™log 3 iV' 



(182) 



We note that all cumulants are of order unity for t ~ log 3 N, which is the sign 
that the distribution of the front position is far from being a trivial Gaussian. 
This makes it particularly interesting. On the other hand, the cumulants are 
proportional to k — t/ log 3 N, which is the sign that the position of the front 
is the result of the sum of k independent random variables, and as such, 
becomes Gaussian when k is very large. The properties of the statistics of the 
front position were investigated in some more details in Ref. [112]. 

Thanks to our discussion in Sec. 2, we see that these results should apply to 
QCD with the relevant substitution of the kernel u and of the parameter N 
according to Tab. 1. 



4-3.2 Numerical simulations 

These results rely on a number of conjectures that no-one has been able to 
prove so far. In order to check our results, let us consider again the model 
introduced in Sec. 2.2.2. The first step to take before being able to apply our 
results to this particular model is to extract from the linear part of Eq. (23) the 
corresponding function cu( 7 ), and then to compute 7c . Setting Ax = At = 1, 
we get 

u>(i) = log (l + A + Pl (e^ - 1) + Pr (e 7 - 1)) , (183) 
and 7c is defined by u;( 7c ) = 7c cj'( 7c ). 

For the purpose of our numerical study, we set 

p l=Pr = 0.1 and A = 0.2. (184) 

Simulated realizations for this set of parameters are shown in Fig. 25. 

From (183), this choice leads to 

7c = 1.352-.. , o/( 7c ) = 0.2553-.., 
u"( lc ) = 0.2267- ••. 

Predictions for all cumulants of the position of the front are obtained by 
replacing the values of these parameters in Eqs. (181), (182). 

Technically, in order to be able to go to very large values of N, we replace the 
full stochastic model by its deterministic mean field approximation u — > (u), 
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Fig. 25. 1000 realizations of the model introduced in Sec. 2.2.2 at two different times 
(dotted lines) , and the average of u over the realizations (full line) . One clearly sees 
that (u) does not keep its shape upon time evolution, which shows that the traveling 
wave property of the FKPP equation is lost due to the stochasticity. 

where the evolution of (u) is given by Eq. (23), in all bins in which the num- 
ber of particles is larger than 10 3 (that is, in the bulk of the front). Whenever 
the number of particles is smaller, we use the full stochastic evolution (21). 
We add an appropriate boundary condition on the interface between the bins 
described by the deterministic equation and the bins described by the stochas- 
tic equation so that the flux of particles is conserved [113]. This setup will be 
called "model I" . Eventually, we shall use the mean field approximation every- 
where except in the rightmost bin (model II): at each time step, a new bin is 
filled immediately on the right of the rightmost nonempty site with a number 
of particles given by a Poisson law of average 9 = N (u(x, t+l)\{u(x, t)}) . We 
checked numerically that this last approximation gives indistinguishable re- 
sults from those obtained within model I as far as the statistics of the position 
of the front is concerned. 

We define the position of the front at time t by 

oo 

X t = Y,u(x,t). (186) 

x=0 

We start at time t — from the initial condition u(x, 0) = 1 for x < 
and u(x, 0) = for x > 0. We evolve it up to time t = log 2 N to get rid 
of subasymptotic effects related to the building up of the asymptotic shape 
of the front, and we measure the mean velocity between times log 2 N and 
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Fig. 26. [From Ref. [Ill]] From top to bottom, the correction to the velocity given 
by the cutoff theory and the cumulants of orders 2 to 5 of the position of the front 
in the stochastic model. The numerical data are compared to our parameter-free 
analytical predictions (181), (182), represented by the dashed line. 

16 x log 2 N. For model I (many stochastic bins), we average the results over 
10 4 such realizations. For model II (only one stochastic bin), we generate 10 5 
such realizations for iV < 10 50 and 10 4 realizations for iV > 10 50 . In all our 
simulations, models I and II give numerically indistinguishable results for the 
values of N where both models were simulated, as can be seen on the figures 
(results for model I are represented by a circle and for model II by a cross). 



Our numerical data for the cumulants is shown in Fig. 26 together with the an- 
alytical predictions obtained from (181), (182) (dotted lines in the figure). We 
see that the numerical simulations get very close to the analytical predictions 
at large N. However, higher-order corrections are presumably still important 
for the lowest values of N displayed in the figure. 

We try to account for these corrections by replacing the factor (log N) / r y c = L 
in the denominator of the expression for the cumulants in Eqs. (181), (182) by 
the ansatz 

31og(logAQ >g(logAQ n „_, 

L cfr = LH Vc + d — — . (187) 

7c log N 

The two first terms in the r.h.s. are suggested by our model. We have added 
two subleading terms which go beyond our theory: a constant term, and a term 
that vanishes at large N. The latter are naturally expected to be among the 
next terms in the asymptotic expansion for large N. We include them in this 
numerical analysis because in the range of N in which we are able to perform 
our numerical simulations, they may still bring a significant contribution. 
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We fit (187) to the numerical data obtained in the framework of model II, 
restricting ourselves to values of N larger than 10 30 . In the fit, each data point 
is weighted by the statistical dispersion of its value in our sample of data. We 
obtain a determination of the values of the free parameters c = —4.26 ± 0.01 
and d — 5.12 ± 0.27, with a good quality of the fit (x 2 /d.o.f ~ 1.15). 

Now we see that with this modification, the results for the cumulants shown 
in figure 26 (full lines), are in excellent agreement with the numerical data 
over the whole range of N. 



5 Application to the computation of QCD scattering amplitudes 

In this section, we shall study the phenomenological relevance of the results 
obtained from the correspondence with statistical physics. There are two as- 
pects that should be discussed. First, we go back to the assumptions that were 
required to go from QCD to reaction-diffusion, and in particular, the hypoth- 
esis of uniformity of the evolution in impact-parameter space. Next, we derive 
new properties of the QCD scattering amplitudes and discuss their impact on 
phenomenology. 

5.1 Relevance of one- dimensional models: impact-parameter correlations 

So far, we have argued that high-energy scattering in QCD at fixed coupling 
and fixed impact parameter is in the universality class of the stochastic FKPP 
equation (Sec. 2), which is an equation with one evolution variable (time 
or rapidity in QCD), and one spatial dimension (x generically, or log A; 2 ~ 
log(l/r 2 ) in QCD). From the very beginning, we have simply discarded the 
impact parameter dependence. It is important to understand that the spatial 
variable and the impact parameter play different roles, and thus, the impact 
parameter may a priori not be accounted for by a two-dimensional extension 
of the FKPP equation. 

There are general arguments to support the assumption that the QCD evolu- 
tion is local enough in impact parameter for the different impact parameters 
to decouple through the rapidity evolution, which we are now going to present. 

Let us start with a single dipole at rest, and bring it gradually to a higher 
rapidity. As was explained in Sec. 2, during this process, this dipole may be 
replaced by two new dipoles, which themselves may split, and so on, eventually 
producing a chain of dipoles. Figure 3 pictures one realization of such a chain. 

According to the splitting rate given in Eq. (1), splittings to smaller-size 
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dipoles are favored, and thus, one expects that the sizes of the dipoles get 
smaller on the average, and that in turn, the successive splittings become more 
local. The dipoles around region "1" and those around region "2" should have 
an independent evolution beyond the stage pictured in Fig. 3: further split- 
tings will not mix in impact parameter space, and thus, the traveling waves 
around these regions should be uncorrelated. For a dipole in region 1 of size 
r to migrate to region 2, it should first split into a dipole whose size is of the 
order of the distance Ab between regions 1 and 2, up to some multiplicative 
factor of order I. (We assume in this discussion that the dipoles in region 2 
relevant to the propagation of the local traveling waves, that is, those which 
are in the bulk of the wave front, also have sizes of order r). Roughly speak- 
ing, the rate of such splittings may be estimated from the dipole splitting 
probability (1): it is of order a(r 2 /(Ab) 2 ) 2 , while the rate of splittings of the 
same dipole into a dipole of similar size in region 1 is of order a. Thus the 
first process is strongly suppressed as soon as regions 1 and 2 are more distant 
than a few units of r. Note that for Ab > 1/Q S , saturation may further reduce 
the emission of the first, large, dipole leading to an even stronger suppression 
of the estimated rate. 

What could also happen is that some larger dipole has, by chance, one of 
its endpoints tuned to the vicinity of the coordinate one is looking at (at 
a distance which is at most |Ar| <C 1/Q S (Y)), and easily produces a large 
number of dipoles there. In this case, the position of the traveling wave at 
that impact parameter would suddenly jump. If such events were frequent 
enough, then they would modify the average wave velocity and thus the one- 
dimensional sFKPP picture. We may give a rough estimate of the rate at which 
dipoles of size smaller than Ar are produced. Assuming local uniformity for 
the distribution n of the emitting dipoles, the rate (per unit of ay) of such 
events can be written 



where we integrate over large dipoles of size r > Ar emitting smaller dipoles 
(of size e < Ar) with a probability d 2 e r 2 /(27re 2 (r — e) 2 ). The factor (e/r ) 2 
accounts for the fact that one endpoint of the dipole of size r has to be in a 
given region of size e in order to emit the dipoles at the right impact parameter. 
To estimate this expression, we first use n(r ) = T(r )/a 2 and approximate T 



The front is replaced by 1 above the saturation scale (for r > 1/Q S ) and by 
an exponentially decaying tail for r < l/Q s . Using r — e ~ r in the emission 
kernel, the integration is then easily performed and one finds a rate whose 
dominant term is 




(188) 



by 






(190) 
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For (Ar) 2 <C (a 2 ) 1 / 72 /Q 2 S , i.e. ahead of the bulk of the front, this term is 
parametrically less than 1 and is in fact of the order of the probability to find 
an object in this region that contributes to the normal evolution of the front 
[111]. Hence there is no extra contribution due to the fact that there are many 
dipoles around at different impact parameters. 

The arguments given here are based on estimates of average numbers of 
dipoles, on typical configurations, and we are not able to account analyti- 
cally for the possible fluctuations. As we have seen through this review, the 
latter often play an important role. As a matter of fact, in the physics of 
disordered systems, rare events sometimes dominate. So before studying the 
phenomenological consequences of the statistical picture of high-energy QCD 
based on a one-dimensional equation, one should check more precisely locality 
of the evolution in impact parameter. 

A numerical check was recently achieved in the case of a toy model that has an 
impact-parameter dependence in Ref. [114]. Let us briefly describe the model. 



5.1.1 A model incorporating an impact-parameter dependence 

In order to arrive at a model that is tractable numerically, we only keep 
one transverse dimension instead of two in 3+1-dimensional QCD. However, 
we cannot consider genuine 2+1-dimensional QCD because we do not wish 
to give up the logarithmic collinear singularities at x 2 = x and x 2 = x\. 
Moreover, QCD with one dimension less has very different properties at high 
energies [115]. Starting from Eq. (1), a splitting rate which complies with our 
requirements is: 



dP 1 


x i\ 




d{ay) 4 


^02! 


\xi 2 \ 



We can further simplify this probability distribution by keeping only its collinear 
and infrared asymptotics (as in Ref. [116]). If \x 02 \ <C |xoi| (or the symmetri- 
cal case |x 12 | <C |#oi|)> the probability reduces to dx 2 /\xQ 2 \ {dx 2 /\xi 2 \ resp.). 
The result of the splitting is a small dipole (x , x 2 ) together with one close in 
size to the parent. So for simplicity we will just add the small dipole to the 
system and leave the parent unchanged. In the infrared region, a dipole of size 
\xq 2 \ ^ l^oi I is emitted with a rate given by the large-|xo2| limit of the above 
probability. The probability laws (1),(191) imply that a second dipole of sim- 
ilar size should be produced while the parent dipole disappears. To retain a 
behavior as close as possible to that in the collinear limit, we will instead just 
generate a single large dipole and maintain the parent. To do this consistently 
one must include a factor of two in the infrared splitting rate, so as not to 
modify the average rate of production of large dipoles. 

In formulating our model precisely, let us focus first on the distribution of the 
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sizes of the participating dipoles. (The simplifying assumptions made above 
enable one to choose the sizes and the impact parameters of the dipoles suc- 
cessively). We call r the modulus of the emitted dipole, r the modulus of its 
parent and Y = ay. The splitting rate (191) reads, in this simplified model 

dP r() ^ r n/ ^r Q dr n/ ^dr 



dY 



= 0(r - r )^V + 6(r - r) — , (192) 



and the original parent dipole is kept. Logarithmic variables are the relevant 
ones here, so we introduce 

p = log 2 (l/r) or r = 2~ p . (193) 

We can thus rewrite the dipole creation rate as 
dP 

-jp = 0(p o -p)2< , -"> \og2 dp + 9(p-p ) log 2 dp. (194) 

To further simplify the model, we discretise the dipole sizes in such a way that 
p is now an integer. This amounts to restricting the dipole sizes to negative 
integer powers of 2. The probability that a dipole at lattice site i (i.e. a dipole 
of size 2~ l ) creates a new dipole at lattice site j is 

dPi^j fpj+i dP n 
dY 



/■«+! d/V^,_ log; ]>i 
L ~7V- 2i < j<i' (195) 



The rates dP i± /dY for a dipole at lattice site % to split to any lattice site j > % 
or j < i respectively are then given by 

rl p. rl P ■ 

3 ~] (196) 
dY ~f^ dY 

where we have restricted the lattice to < % < L, for obvious reasons related 
to the numerical implementation. 

Now we have to address the question of the impact parameter of the emitted 
dipole. In QCD, the collinear dipoles are produced near the endpoints of the 
parent dipoles. Let us take a parent of size r at impact parameter b . We set 
the emitted dipole (size r) at the impact parameter b such that 

where s has uniform probability between and 1. It is introduced to obtain 
a continuous distribution of the impact parameter unaffected by the discreti- 
sation of r. This prescription is quite arbitrary in its details, but the latter 
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do not influence significantly the physical observables. Each of the two signs 
that appear in the above expression is chosen to be either + or — with equal 
weights. We apply the same prescription when the emitted dipole is larger 
than its parent. 

Scattering amplitude 

We have explained above (see Sec. 2.1) that in QCD, the scattering amplitude 
of an elementary probe dipole of size r\ = 2~ % with a dipole in an evolved 
Fock state is proportional to the number of objects which have a size of the 
same order of magnitude and which sit in a region of size of order around 
the impact point of the probe dipole. Since in our case, the sizes are discrete, 
the amplitude is just given, up to a factor, by the number of dipoles that are 
exactly in the same bin of size as the probe, namely 

T(i, b ) =a s 2 x #{dipoles of size 2~ l 

at impact parameter b satisfying \b — b \ < rj/2}. (198) 

Saturation 

We now have to enforce unitarity, that is the condition 

T{i,b)<\ (199) 

for any i and b. This condition is expected to hold due to gluon saturation in 
QCD. However, saturation is not included in the original dipole model. The 
simplest choice is to veto splittings that would locally drive the amplitude to 
values larger than 1. In practice, for each splitting that gives birth to a new 
dipole of size % at impact parameter b, we compute T(i, b) and T(i, b ± n/2), 
and throw away the produced dipole whenever one of these numbers gets larger 
than one. 

Given the definition of the amplitude T, this saturation rule implies that 
there is a maximum number of objects in each bin of size and at each impact 
parameter, which is equal to iV sa t = l/etg- 

5.1.2 Numerical results 

We take as an initial condition a number iV sat of dipoles of size 1 (i — 0), 
uniformly distributed in impact parameter between — r /2 and r /2. The im- 
pact parameters bj that are considered are respectively 0, 10 -6 , 10 -4 , 10~ 2 
and 10 _1 . The number of events generated is typically 10 4 , which allows one 
to measure the mean and variance of the position of the traveling waves to a 
sufficient accuracy. 
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We have checked that at each impact parameter, we get traveling waves whose 
positions grow linearly with rapidity at a velocity less than the expected mean- 
field velocity for this model. iV sat was varied from 10 to 200. 

Fig. 27 represents the correlations between the positions of the wave fronts at 
different impact parameters in the AIP model, defined as 



{p s (YM)Ps{YM)) - ( Ps (YM)){Ps(YM)}- (200) 



We set N sat to 25 in that figure, but we also repeated the analysis for different 
values of iV sa t between 10 and 200. 

We see very clearly the successive decouplings of the different impact parame- 
ters in Fig. 27, from the most distant to the closest one, as rapidity increases. 
Indeed, the correlation functions flatten after some given rapidity depending 
on the difference in the probed impact parameters, which means that the evo- 
lutions decouple. This decoupling is expected as soon as the traveling wave 
front reaches dipole sizes which are smaller than the distance between the 
probed impact parameters, i.e. at Y such that \b 2 — b±\ ~ 1/Q S (Y) = 2~ ps ( Y \ 
From the data for p 8 (Y), we can estimate quantitatively the values of the ra- 
pidities at which the traveling waves decouple between the different impact 
parameters. (It is enough to invert the above formula for the relevant values of 
62 — bi). These rapidities are denoted by a cross in Fig. 27 for the considered 
impact parameter differences. Our numerical results for the correlations are 
nicely consistent with this estimate, since the correlations start to saturate to 
a constant value precisely on the right of each such cross. 

We conclude that the different impact parameters indeed decouple, as was 
expected from a naive analytical estimate. What is true for our toy model 
should go over to full QCD, since we have included the main features of QCD. 
When looking at the numerical data more carefully however, it turns out that 
the model with impact parameter does not reduce exactly to a supposedly 
equivalent one-dimensional model of the sFKPP type. This is a point that 
would deserve more work. We refer the reader to Ref. [114] for all details of 
our numerical investigations. 

An attempt to build a complete picture of high-energy QCD that includes the 
impact parameter was made in Ref. [117], but it relies on some more conjec- 
tures, that are difficult to prove. Finding a mathematically sound formulation 
remains a challenge. 
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Fig. 27. Correlations of the positions of the traveling wave fronts between different 
impact parameters in the toy model of Sec. 5.1. The points where the correlations 
flatten correspond to the decoupling of the waves in the corresponding regions of 
impact parameter. 



5.2 Traveling waves, geometric scaling, and consequences of the noise 



As was stated in the Introduction, the initially unplanned opportunity to 
collect data in the high-energy regime of deep-inelastic scattering at HERA 
triggered a renewed interest in small-x physics among phenomenologists. The 
major discoveries in this regime is the (unexpected) important fraction of 
diffractive events, and a new scaling, geometric scaling, featured by total (and 
even semi-inclusive) cross sections (see Fig. 1). 

In order to deal theoretically with the small- a: regime, one needs new fac- 
torization theorems in order to single out the elements of the cross sections 
that are computable in perturbation theory. High-energy, or /cj_-factorization, 
[118,119,120] is the appropriate tool. A practical way to implement /cj_-factorization 
is the color dipole model presented in Sec. 2. 
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5.2.1 Dipole models and geometric scaling 



The main observable measured at HERA is the proton structure function F 2 . 
It is proportional to the sum of the virtual photon-proton cross section for a 
transversely and longitudinally polarized photon respectively. 

A bare photon has no hadronic interactions, since it does not carry any color 
charge. However, it may easily fluctuate into a quark- ant iquark pair, overall 
color-neutral, thus forming a color dipole. Subsequently, these dipoles will 
interact with the target proton. This picture is represented by the following 
equations: 

^ = i£;(-+-). (201) 

°~t,l{ x iQ 2 ) = j dzd 2 r\^ T:L (z,r,Q 2 )\ 2 a dipo i c (x,r). 

Here, o~t,l are the photon-proton cross sections for transversly and longitu- 
dinally polarized virtual photons. \&t,l are light-cone wavefunctions for 7*, 
computable within QED (see, e.g., Ref. [27] for explicit expressions to lowest 
order in a em ). Furthermore, o"di po ie(^, t) is the cross-section for dipole-proton 
scattering (for a dipole of transverse size r), and encodes all the information 
about hadronic interactions (including unitarization effects). This cross sec- 
tion is related to the amplitude A discussed so far by an integration over the 
impact parameter. (Actually, A was the forward elastic amplitude; the optical 
theorem relates it to the total cross section). 

In Ref. [27,28], the dipole cross-section was modeled as 

^dipoic(^, r) = (To (l - e^^ 4 ) , (202) 

where a is a hadronic cross-section: It stems from the integration over the 
impact parameter, when the impact parameter dependence is supposed to be 
uniform over a disk of radius ~ \fo\)- Qs{%) plays the role of the saturation 
momentum, parametrized as Q 2 (x) = (x /x) x x 1 GeV 2 . Note that, by con- 
struction, this cross section only depends on the combined variable r 2 Q 2 (x) 
instead of r and x separately. This property is transmitted to the measured 
photon cross sections <Tt,l(x, Q 2 ), which then depend on Q 2 /Q 2 (x) only (this 
scaling is slightly violated by the masses of the quarks). This is geometric 
scaling, predicted to be a feature of the solutions to the BK equation at large 
rapidity. 

Historically, geometric scaling was discovered first in the data (see Ref. [29]), 
after Golec-Biernat and Wiisthoff (GBW) had written down their model: The 
latter happened to feature this scaling (up to small violations induced by 
the quark masses). There was no apparent need for finite rapidity scaling 
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violations in the first HERA data. However, later analysis revealed that a 
significant amount of explicit scaling violations in the dipole cross section, 
predicted by the BK equation, were actually required by more accurate data. 

A now popular model that describes the HERA data in a way that takes a 
better account of the subasymptotics, beyond the GBW model, was formu- 
lated in Ref. [121]. The dipole scattering cross section reads <7di P o\e(x, r) = 
27rR 2 J\f(y,rQ s ), with 



where Q s = Q s (x) = (x /x) x ^ 2 GeV. The expression for the cross section 
for r small compared to 2/Q s corresponds to the solution of the BK equa- 
tion (compare to Eq. (143) with the help of Tab. 1), in which we substituted 
^(7c) = x(lc) an d u/'(7 c ) = x"(lc) by the parameters A and n that we subse- 
quently fit to the data. The expression in the second line also has the correct 
functional form for r 3> 2/Q s , as obtained by solving the BK equation [32]. 
This is strictly valid only to leading-order accuracy, but here it is used merely 
as a convenient interpolation towards the 'black disk' limit M — 1. (The details 
of this interpolation are unimportant for the calculation of <7 7 * p .) The coef- 
ficients a and b are determined uniquely from the condition that J\f(rQ s , Y) 
and its slope be continuous at rQ s = 2. The overall factor Af in the first line 
of Eq. (203) is ambiguous, reflecting an ambiguity in the definition of Q s . This 
model fits well all HERA data for structure functions, in the range x < 10 -2 . 
All details may be found in Ref. [121]. 

The model explicitely breaks geometric scaling. However, effectively, geometric 
scaling remains a fairly good symmetry of the model, as required by the data. 
The small finite-rapidity scaling violations are needed to describe accurately 
the high-precision HERA data. 

The model may also accomodate less inclusive observables, such as diffraction 
[122]. It has been improved recently by including heavy quarks [123] (The cru- 
cial need for taking account of the charm quark was emphasized in Ref. [124]). 
An impact-parameter dependence was also introduced [125,126,127] that was 
already missing in the GBW model. 

The range of validity of dipole models has been re-examined recently [128]. 
5.2.2 Diffusive scaling 

At still higher energies, according to the discussion of Sec. 4, one expects the 
saturation scale to acquire a dispersion from event to event that scales with 




for rQ s <2, 
for rQ s > 2, 



(203) 
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the rapidity like y/ay when rapidity increases. Although this dispersion is not 
an observable since there is no way to measure the saturation scale of an 
individual event, it manifests itself in the total cross section in the form of a 
new scaling, different from geometric scaling. 

The physical amplitude for the scattering of a dipole of size r off some target 
is given by the average of all realizations of the evolution at a given y: 

A(y,r) = (T(r))\ y . (204) 

For large enough rapidities and small enough a s , these realizations are expo- 
nentially decaying fronts in the variable p = log(l/r 2 ), fully characterized by 
a stochastic saturation scale, or rather its logarithm p s = log Q 2 s (y). For the 
purpose of the present discussion, it may be approximated in the same way 
as in Eq. (189), namely 

T(p) = 9( Ps - p) + 9{p - p s )e-^ p - p °\ (205) 

The statistics of p s is given by Eqs. (181), (182) (up to the replacements sug- 
gested in Tab. 1 to go from a generic reaction-diffusion to QCD). At ultrahigh 
energies (and very small a s ), it is essentially a Gaussian centered at 

(Ps) = ^' {lc) ,) ay (206) 

V 7e 2(log(l/a2) + 31oglog(l/^)) 2 y 



and of variance 



° 2 = (Pi) - (Ps) 2 = %\, 2 ^ y. (207) 
3 log [l/al) 

The scattering amplitude may be expressed by the simple formula 



MV,P) = ^J dpsT(p)\yew Jf ))2 ) • (208) 

The most remarkable feature of this amplitude is the scaling form for A that 
it yields: 

A(v,p) = a( ' -<"•(»)> V (209) 

VvWto?(iK)/ 

This equation may be obtained by performing the integration in Eq. (208) 
after the replacement of T by its approximation (205). This scaling obviously 
violates geometric scaling: If the latter scaling were satisfied, then A would be 
a function of p — (p s (y)) only. 

Mueller and Shoshi had already noted that geometric scaling had to be vio- 
lated beyond the BK equation in Ref. [39]. However, the square root in the 
denominator of the scaling variable in Eq. (209) was missing because their 
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approach was relying on mean field throughout, thus missing the stochastic 
nature of the evolution. 



This new scaling is a firm prediction of the correspondence with statistical 
physics. However, it may not be tested at particle colliders in a simple way. 
Let us work out the order of magnitude of the rapidity needed for the differ- 
ent effects (saturation, geometric scaling, diffusive scaling) to show up. The 
rapidity that is needed to reach saturation is roughly 

2/BFKL _ TTn ■ (210) 

The BK picture is expected to be valid until the asymptotic exponential shape 
of the front has diffused down to the point where the amplitude becomes of the 
order of o? s . This additional rapidity needed to get to the regime of geometric 
scaling is thus given by Eq. (157) once the appropriate replacements have been 
done 



1 

2/BK 



log(l/«, 2 ) 



2«x"(7c) L 7c J ' (211) 
and finally, the effect of the fluctuations of the saturation scale gets important 
at the rapidity 

3iog 3 (i/«;) ,„,„, 

ct ~ carVM ' <212) 
The relevant parameters read, in QCD, 

7c = 0.627549, X {lc) = 3.0645, X " (7c) = 48.5176. (213) 

For some realistic strong coupling constant, a s ~ 0.2, we get 

t/BFKL ~ 6.07879 , y BK ~ 1.41965 , y Ruct ~ 0.348244. (214) 

Given that rapidities in the small- x regime at HERA were of the order of 
10, and will be of the order of 15 at the LHC, these figures may give us the 
hope that we may observe these effects. However, the values of the rapidity 
that delimitates the different regimes are largely underestimated given that 
they rely on the leading-order BFKL kernel, which predicts a much too large 
growth of the cross section with the rapidity and a too fast diffusion (see the 
large value of x"(lc))- One also has to keep in mind that the former estimates 
should only hold for very small values of a s . 

Furthermore, the effect of the running coupling, which should be taken into 
account in any detailed phenomenological study, is expected to still reduce the 
effects of the fluctuations [129]. 

Nevertheless, the effect of diffusive scaling (i.e. of the event-by-event fluctua- 
tions of the saturation scale) on observables has already been investigated in 



82 



some detail by several groups. Diffractive amplitudes were studied in Ref. [130]. 
The ratio of the gluon distribution in a nucleus to the same quantity in a pro- 
ton was computed in Ref. [131]. 



6 Conclusion and outlook 

We have reviewed a peculiar way of viewing high-energy scattering in QCD, 
based on the physics of the parton model, and its strong similarities with 
reaction-diffusion processes (Sec. 2). The correspondence is best summarized 
in the mapping of Tab. 1. We have seen that the equations that describe 
the dynamics of these processes are in the universality class of the stochastic 
FKPP equation, and admit traveling-wave solutions whose features are likely 
to be universal, in such a way that a study of simple reaction-diffusion-like 
models may lead to exact asymptotic results also for QCD scattering ampli- 
tudes. Understanding the very mechanism of traveling wave formation and 
front propagation was crucial to see how the universality may come about 
(see Sec. 4). 

In zero-dimensional stochastic models, we could perform exact calculations 
and get analytical results within different formulations (Sec. 3). We understood 
that analyzing the structure of single events was technically much simpler if 
one wants to get leading orders at large N (= l/a^), since in individual realiza- 
tions, one may factorize the fluctuating part from the nonlinear effects. Thanks 
to this observation, in one-dimensional models which admit realizations in the 
form of stochastic traveling waves, we could also get precise analytical results 
on the form and shape of the traveling waves, which are presumably exact 
asymptotically (Sec. 4). Universality enables one to make statements on the 
form of the QCD scattering amplitudes at very high energies. These state- 
ments turn into firm phenomenological predictions (Sec. 5), which however do 
not seem to be testable at colliders in the near future. Nevertheless, getting 
new analytical results for QCD in some limit is always an interesting achieve- 
ment, given the complexity of the theory. Furthermore, while our analytical 
results only apply for exponentially small a s (log(l/a^) 1), the picture 
itself should be valid in the whole perturbative range, namely for C 1. 

There are still many open questions. On the statistical physics side, the statis- 
tics of the front position that we have found has not been derived rigorously, 
but rather guessed, and rely on many quite ad hoc conjectures. We got confi- 
dence on the validity of our conjectures on the basis of numerical simulations. 
Moreover, although we expect universality up to corrections of order 1/N (that 
is to say 0(a 2 s ) in QCD), we could only get analytical expressions relative to 
the cumulants of the position of the front for the first terms in an expansion 
in powers of 1/ log N, which requires much larger values of N to be valid. On 
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a more general footing, the sFKPP equation seems to describe many physi- 
cal, chemical or biological problems (in particular population evolution with 
selection in evolutionary biology). We have also found recently an explicit 
analogy with the theory of spin glasses [132]. This large universality is maybe 
the strongest incentive to try and find more accurate solutions to that kind of 
equations. 

On the QCD side, the correspondence with reaction-diffusion processes strongly 
relies on the assumption that there is saturation of some form of the quark 
and gluon densities in the hadronic wave functions. While this is a reasonable 
guess that few experts would challenge, it is clear that we cannot consider 
that the problem is solved before the saturation mechanism at work in QCD 
has been exhibited. QCD is formulated as a quantum field theory. To see the 
similarity with reaction-diffusion, we basically needed to translate it into the 
parton model first. It would be better to recover the results of Sec. 4 (and 
hopefully get more) directly from field theory [133], as one could do it in 
the zero-dimensional model introduced in Sec. 3. This requires to understand 
the strong field regime of field theory. This is an exciting challenge for both 
particle physicists and statistical physicists. 
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